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ABSTRACT 


A  shell  formulation  was  developed  from  a  three-dimensional  solid.  The  shell 
element  has  four  comer  nodes  at  which  there  are  three  displacements  and  three  rotations  as 
nodal  degrees  of  freedom,  and  includes  both  transverse  shear  and  transverse  normal 
deformations.  The  element  utilizes  reduced  integration  along  the  in-plane  axes  and  full 
integration  along  the  transverse  axis.  The  formulation  incorporates  the  Gurson  constitutive 
model  for  void  growth  and  plastic  deformation.  An  algorithm  for  stable  solutions  of  the 
nonlinear  constitutive  equations  is  also  developed.  Hourglass  mode  control  is  provided  by 
adding  a  small  fraction  of  internal  force  determined  through  full  integration  along  the  in¬ 
plane  axes  and  reduced  integration  along  the  transverse  axis.  Implementation  into  a  research 
finite  element  program  is  discussed.  Numerical  examples  are  provided  to  verify  the  accuracy 
of  the  new  element  and  to  show  the  importance  of  the  transverse  normal  stress,  void  effects 
on  plastic  strain,  and  the  necessity  of  applying  a  drilling  moment 
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I.  INTRODUCTION 


Since  shell  structures  are  efficient  load-carrying  structural  members,  they  have  been 
dominant  in  most  structural  applications.  However,  the  finite  element  formulation  of  shells 
poses  some  difficulty.  As  a  result,  extensive  study  has  been  devoted  to  developing  better 
shell  elements.  Because  of  the  abundance  of  papers  on  this  subject,  no  attempt  is  made  here 
to  summarize  all  of  them.  Some  of  the  related  past  work  is  given  in  references  [1-12]. 

Void  growth  and  nucleation  can  have  a  significant  effect  on  plastic  flow  [13].  Since 
voids  act  as  stress  concentrators,  the  overall  effect  is  to  reduce  the  stress  under  plastic  flow, 
and  increase  the  plastic  strain  [14-16].  The  model  proposed  by  Gurson  [17]  appears  to  have 
been  adopted  as  the  standard  for  incorporating  void  growth  and  nucleation  effects  into  a 
numerical  solution  of  elasto-plastic  problems.  Previous  work  has  been  done  on  improving 
the  efficiency  and  accuracy  of  plasticity  computations,  applying  plasticity  to  plate/shell 
elements,  and  incorporating  void  growth  and  nucleation  effects  in  solid  elements  [18-20]. 
The  element  presented  in  this  paper  incorporates  a  modified  version  of  the  algorithm  for 
three-dimensional  solids  proposed  by  Aravas  [21]  into  a  shell  element.  This  shell  element 
assumes  a  modified  plane-stress  condition,  and  utilizes  both  the  hydrostatic  pressure  and  the 
deviatric  stress.  Due  to  the  importance  of  the  hydrostatic  pressure  on  void  growth  and 
nucleation,  this  element  includes  the  transverse  normal  stress. 

The  algorithm  proposed  by  Aravas  is  not  unconditionally  stable  when  applied  to  this 
modified  plane  stress  condition  [21].  This  paper  also  presents  modifications  to  the  original 
algorithm,  which  greatly  enhance  the  stability  of  the  solution,  at  the  cost  of  a  slight  decrease 
in  computational  efficiency.  This  algorithm  also  allows  for  more  complicated  work- 
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hardening  profiles  than  the  typical  exponent  law  (cr  =  Ken).  Full  modeling  of  the  entire 
stress-strain  curve  is  accomplished  by  a  piece-wise  linear  approximation.  While  this  method 
is  not  particularly  useful  for  analytic  approaches,  it  appears  to  be  helpful  in  strictly  numerical 
solutions.  The  stress-strain  relationship  may  be  taken  directly  off  the  results  of  a  standard 
tensile  test. 

For  thick  shell  applications,  the  transverse  normal  stress  and  strain  cannot  be 
neglected.  The  transverse  normal  stress  and  strain  affect  both  the  hydrostatic  pressure  and 
the  deviatoric  stress,  which  in  turn  affect  plastic  deformation  and  void  growth  and 
nucleation.  The  element  presented  here  extends  the  typical  use  of  the  drilling  degree  of 
freedom  (DOF),  as  outlined  in  Hughes  and  Brezzi  [22],  to  incorporate  the  effects  of 
transverse  normal  strain.  The  drilling  DOF  is  important  in  transmitting  stress  and  strain  to 
elements  across  sharp  bends  and  curves,  and  is  therefore  already  included  in  a  variety  of  shell 
elements  [2,4,6,7,10].  In  addition  to  this  traditional  usage,  the  present  element  utilizes  the 
drilling  DOF  to  compute  the  transverse  normal  deformation.  The  importance  of  including 
the  transverse  normal  deformation  is  outlined  in  Essenburg  [23]. 

The  present  study  formulates  a  shell  element  for  transient  analysis.  The  element  can 
have  elasto-plastic  deformation  with  void  growth  and  nucleation.  Gurson's  void  model  is 
used  as  a  basis  for  the  void  constitutive  model.  The  shell  element  includes  both  transverse 
shear  deformation  and  the  transverse  normal  deformation  for  thick  shell  applications.  The 
drilling  degree  of  freedom  is  used  for  computing  the  deformation  through  the  thickness  of 
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a  thick  shell.  An  algorithm  for  stable  solutions  of  the  nonlinear  constitutive  equations  is  also 
developed. 

Some  example  problems  are  presented  to  evaluate  the  formulation  and  to  investigate 
the  effects  of  the  transverse  normal  strain  for  thick  shells  in  association  with  elasto-plastic 
deformation,  including  void  effects.  The  sections  that  follow  are:  the  formulation  of  the  shell 
element,  discussion  of  the  elasto-plastic  constitutive  equation  including  void  nucleation  and 
growth,  the  stable  solution  algorithm  including  spurious  (hourglass)  mode  control,  numerical 
examples  and  discussion  of  results,  and  conclusions. 
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H.  FINITE  ELEMENT  FORMULATION 

A.  GEOMETRY 

A  point  in  a  shell  structure  can  be  expressed  by  a  vector  sum  of  two  vectors.  The  first 
vector  is  a  position  vector  from  the  origin  of  the  global  coordinate  system  to  a  point  on  a 
reference  surface  of  the  shell  element.  The  second  vector  is  a  position  vector  from  this 
reference  surface  to  the  point  under  consideration.  The  surface  that  spans  the  center  of  the 
transverse  axis  is  used  as  the  reference  surface  in  this  formulation,  although  any  surface 
would  suffice.  The  first  vector  terminates  at  the  reference  surface  directly  below  the  point 
in  question.  The  second  vector  is  then  the  normal  from  the  reference  surface  that  intersects 
the  desired  point  Figure  1  shows  this  relationship.  Two  shape  functions  are  used  to  describe 


Figure  1.  Element  Cross  Section. 
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a  position  in  the  element;  is  the  two  dimensional  shape  function  in  the  plane,  and 
H  k  is  the  one  dimensional  shape  function  along  the  C,  axis,  where  (^,r|,Q  describes  a  point 
in  the  natural  coordinate  system.  A  generic  point  in  the  shell  may  now  be  described  in  terms 
of  the  position  vectors  of  the  nodes  and  the  shape  functions: 

x,(4.’>,0  =  'ZNt(4.’J)  X1;  +tlNla,n)Hk(()Vks,  o- 1.2.3)  (1) 

k=l  k= 1 


where  xl  is  the  position  vector  of  node  A:  in  the  reference  surface;  V%  is  the  unit  vector  at 
the  node  k,  and  n  is  the  number  of  nodes  per  element.  In  the  present  formulation,  a  four-node 
shell  element  is  considered.  The  unit  vector  Vk3,  is  defined  as: 

k  tfr-ti)*- 

3i  ||(*?)<<,p  1  j 


where  top  and  bottom  indicate  the  top  and  bottom  surfaces  of  the  shell,  and  ||  ||  denotes  the 
Euclidean  norm.  The  one-dimensional  shape  function  H*  is  expressed  as: 


Hk(0  = 


±(i+<rxi-o-;j(i-<rxi+o 


bottom 


(3) 


in  which  £  indicates  the  location  of  the  reference  surface  and  varied  from  -1  to  1(^  =  0 
denotes  the  mid-surface).  The  two-dimensional  shape  function  is  expressed  as: 

N2=i(l+S)(l-V) 

N3  =  i(l+Z)(l  +  V) 

N<  =  i(l-t)(l+Tj) 
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B.  DISPLACEMENT 


The  displacement  field  in  a  shell  can  be  written  as: 
mtt.rj.O  =  I>Y^  uf  +  XN^^T?) Hk(C)(-vW  +VlJi0l2+Vk3lfi)  0  =  1, 2, 3)  (5) 

k=\  k=l 

in  which  u\  is  the  displacement  along  the  x\  axis,  ul  is  the  nodal  displacement  at  the  node 
k,  and  unit  vectors  Vn  and  ]/k2i  lie  along  the  reference  surface.  Vu ,  Vi  and  yk3i  are  mutually 
perpendicular.  0k ,  0k  and  0k3  are  rotational  degrees  of  freedom  along  the  unit  vectors 

Vu ,  Vv  and  ykf ,  respectively.  The  right-hand  rule  is  assumed  for  the  positive  direction  of 
each  rotation.  Figure  2  illustrates  the  relationship  among  these  vectors. 


Figure  2.  Displacement  Vector  Orientation. 
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0\  and  0k2  are  the  bending  rotations,  while  0k3  is  the  drilling  rotational  degree  of 

freedom.  The  role  of  Qk3  can  be  clarified  by  considering  a  flat  plate  parallel  to  the  xrx2 
plane.  Now  Eq.  (5)  can  be  rewritten  for  the  transverse  displacement  as: 

«3  =  J  Nkuk  +  Y  NkHk0k  (6) 

*=1  i=l 

Equation  (6)  demonstrates  that  the  transverse  deformation  varies  through  the  plate  thickness 
(i.e.  along  the  C,  axis)  with  6k .  In  this  way,  the  transverse  normal  deformation  is  included 
in  this  formulation,  along  with  the  transverse  shear  deformations. 

C.  COORDINATE  TRANSFORMATION 

Combining  the  three  unit  direction  vectors  into  matrix  [Tp]  provides  the  rotation 
transformation  matrix,  or  matrix  of  direction  cosines,  as  shown  in  Eq.  (7).  [Tp]'1  is  used  to 
transform  the  nodal  degrees  of  freedom  from  the  global  coordinate  system  to  local 
coordinates,  as  shown  in  Eq.  (8),  where  k  is  the  node  number.  The  components  of  {d}  at  the 
four  nodes  of  an  element  are  shown  in  Eq.  (9).  Once  the  internal  force  vector  is  generated 
in  local  coordinates,  [Tp]  is  used  to  transform  the  local  vectors  back  into  global  coordinates. 
This  procedure  will  be  discussed  later. 

[Tp  ]~|  Vi  V2  V3  1(3x3)  (7) 
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{d}-{d'd2d3d4}T  (9) 

The  strain  transformation  matrix,  [T],  is  used  to  transform  calculated  strain  from  the 
global  coordinate  system  to  local  coordinates.  Transforming  the  resulting  stress  from  local 
coordinates  to  the  global  coordinate  system  would  normally  require  using  [T]'1,  but  since  [T] 
is  orthogonal,  [T]  **  =  [T] T,  where  [T] T  is  the  transpose  of  [T].  This  property  negates  the 
requirement  to  invert  a  six  by  six  matrix.  For  a  detailed  derivation  of  these  transformations, 
refer  to  Cook  [24].  The  strain  transformation  matrix  is  explicitly  defined  in  Eq.  (10),  where 
Vjj  is  the  cosine  of  the  direction  vector  V;  in  the  xj  direction. 
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D.  STRAIN  DISPLACEMENT  RELATION 


The  six  components  of  the  strain  tensor  are  computed  from  Eq.  (5)  by  taking  its 
derivative  with  respect  to  the  xj  axis.  In  matrix  form,  the  result  for  a  four-node  element  is: 


w-ww 


W-{ 


£11  £22  £33  7 12  7  23  T 1 3 


where 


[b]=[[b'][b!][b3][b4]] 


The  detailed  expression  for  [Bk]  is: 


[b*]- 


~dNk 

dxj 

0 

0 

-§^2, 

gjVn 

gkVk3I 

0 

dNk 

dX2 

0 

-g\Vk22 

gXVn 

g\Vk32 

0 

0 

dNk 

8x3 

~  g\  V23 

g\VkU 

g\V\3 

ML 

dx2 

Mi 

dxi 

0 

-gk2Vk2,-g)Vk22 

gk2Vn  + gk,Vkn 

gk2Vk31  +  gjVk32 

0 

dNk 

dxs 

SNk 

8X2 

-glV^-glV^ 

gkkVh  +  gkV^3 

glV^  +  glV^ 

dNk 

_dx3 

0 

ml 

3xi 

-g)Vk21-gk,Vk23 

gkVkU+gkjVk,3 

gk3Vk3,  +  g)Vk33_ 

in  which 


(ID 

(12) 


(13) 


(14) 


*  d  Hk 


dXi 


dxt 


The  vector  {d^  }  is  defined  as: 


{dkh{uk,U2Uk3e1®  2®3  } 


(15) 


(16) 
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where  u\  is  the  displacement  along  the  xj  direction  at  node  k,  and  ©f  is  the  rotational 
displacement  about  the  global  Xj  axis  at  node  k.  The  matrix  [B1']  must  be  calculated  for  each 
integration  point. 


E.  JACOBIAN  MATRIX 


Computing  the  derivatives  and  requires  the  Jacobian  matrix,  defined  as 


[j]= 


XI, f  X24  X3.Z 

X 1 ,7]  X2/1  X3jj 

l  Xl,C  X2{  X34 


(17) 


where 


r  -^Xi_spdN  ,  y  nkv'k 

X‘4  2-t  5;  Xi  +  2-1  3£  H  V#  (l 

k=i  og  kmI  eg 

(18) 

-  | +  0  =  1.2, 3) 

V*]  k=J  U7]  k=1  07] 

<19) 

r  -  ^ Xi  -  V  \Tk  ^Hk  r/k  /j-1  ?  3) 

v"ar'Ir  as  v“ (  X  ) 

(20) 

[R]  is  defined  as  the  inverse  of  the  Jacobian  matrix,  [J]  \  Then  the  required  partial 

derivatives  are  defined  as: 

(21) 

(22) 
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F.  STRESS-STRAIN  RELATIONSHIP 

The  strain  calculated  in  Eq.  (1 1)  is  in  the  global  coordinate,  and  is  transformed  to  a 
local  coordinate  system  using 

{  £ local  }  [T]{£eloU}  (23) 

where  [7]  is  defined  in  Eq.  (10).  Stress  is  calculated  from  the  strain  in  the  local  coordinate 
system  using  the  plane-strain  assumption  for  the  in-plane  stress  components: 

E  (  ,  \  Txy  G  Y  x 


Gx  7  2 

1-V 

E 


(£x  +  V£y) 


xy 


Gy  2(v&x~^£y) 

1-V 


Tyz  =  K  G  y 


yz 


(24) 


Gz  zh  Sz 


Txz-K-G  yx 


where  K  is  the  shear  correction  factor,  E  is  the  elastic  modulus,  G  is  the  shear  modulus,  and 
v  is  Poisson’s  ratio.  The  resulting  stresses  are  in  the  local  coordinate  system  and  are 
converted  to  the  global  coordinate  system  with 

W=[T]T{^}  (25) 

{^}-{a,a,a,r„r„r„l  (26) 


where  [T]  is  from  Eq.  (10). 


DRILLING  MOMENTS 


Let  Mf  be  the  transverse  deflection  as  defined  in  Eq.  (6).  The  work  done  by  a 
pressure  loading,  p,  on  an  element  is: 

W=[uiPdA  (27) 


12 


where  Ae  is  the  element  area.  Substituting  Eq.  (6)  into  Eq.  (27),  the  work  is  now: 

W  =  {  X,  f  l  [N]p  dA  +  {  0,  f  l  [n‘H‘  ]p  dA  (28) 

The  first  term  gives  the  conventional  forces  at  each  node,  while  the  second  term  yields  the 
new  nodal  load,  which  will  be  called  the  Drilling  Moment  (DM).  For  example,  if  p  is  a 
concentrated  force  P  at  node  n,  then  the  drilling  moment  associated  with  6"  becomes 
'ACtP,  where  t  is  the  shell  thickness  and  the  mid-plane  is  the  reference  plane. 

When  P  is  applied  at  the  top  plane,  £  equals  one,  and  the  drilling  moment  is  ]/2tP . 
If  P  is  applied  at  the  bottom  plane  (still  with  a  positive  loading  direction),  the  drilling 
moment  is  -]/2tP .  That  is,  the  load  results  in  compression  or  tension  in  the  transverse 
normal  stress  depending  on  whether  the  loading  is  applied  to  the  top  or  bottom  surface  of  the 
shell.  The  transverse  normal  stress  is  assumed  constant  through  the  shell  thickness  for  elastic 
deformation.  However,  as  plastic  deformation  progresses,  the  transverse  normal  stress 
becomes  non-uniform  through  the  thickness  in  order  to  satisfy  the  yield  function. 

The  drilling  moment  direction  coincides  with  the  direction  of  drilling  rotation,  but 
affects  the  transverse  normal  stress.  The  drilling  moment  is  used  as  a  mathematical 
convenience,  and  thus,  lacks  a  physical  interpretation. 
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H.  INTERNAL  FORCE,  MASS  AND  THEIR  ASSEMBLY 


The  stress  resulting  from  Eq.  (25)  is  then  converted  to  an  internal  force  vector  and 
summed  over  all  integration  points  using 

111  nx  ny  n z 

{ fint}=\  J  J  W7  {<*  =  Z  2  [B]T  {cr}wiWjWk|  J|  (29) 

-/  -/  -/  i=l  j=l  k=l 

where  |  J  |  is  the  determinant  of  the  Jacobian  matrix,  nx,  ny,  and  nz  are  the  number  of 

integration  points  in  the  £  7,  and  C,  directions,  respectively,  and  wj,  wj,  and  wfc,  are  the  Gauss 

weights  related  to  those  integration  points.  The  resulting  force  vector,  {fiM  },  must  have  its 

components  transformed  back  to  the  global  coordinate  system  using: 

1  0  0  0  0  0 

0  1  0  0  0  0 

0  0  1  0  0  0. 

0  0  0  (30) 

0  0  0  [tJ 

0  0  0 

A  lumped  mass  method  is  used  for  the  element  mass  matrix.  The  matrix  for  each 
element  is  diagonal,  with  equal  diagonal  elements: 

me  0  0  ...  0 

0  nu  9  ...  0 

0  0  nu  9 

0  0  0  ...  me 

where 
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with  n  =  4  in  this  four-node  element.  Using  a  diagonal  mass  matrix  greatly  simplifies  time 
integration,  since  inverting  it  is  a  trivial  task  requiring  little  computation  time.  The  internal 
force  vector  and  element  mass  vector  (each  24  by  1)  are  then  assembled  into  the 
corresponding  system  vectors. 

I.  EXPLICIT  TIME  INTEGRATION 

The  use  of  internal  force  vectors  and  explicit  time  integration  negates  the  need  to 
explicitly  form  the  system  stiffness  matrices.  The  acceleration  vector  is  computed  from 

(33) 

where  {u  }  is  the  system  acceleration  vector,  [M]  is  the  system  mass  matrix,  { F ext }  is  the 
system  external  force  vector,  and  superscript  t  denotes  the  time  step.  Of  course,  the  system 
mass  matrix  in  Eq.  (33)  is  simply  symbolic,  since  a  system  mass  vector,  formed  from  the 
diagonal  of  the  mass  matrix,  is  used  in  actual  computation.  Velocity  and  displacement  are 
then  found  using 

{up  =  {u}"‘t  +  {£?}'  At 
{u}‘+*  =  {u}‘  +  {c/}'+f  At 
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HI.  DAMAGE  CONSTITUTIVE  EQUATIONS 
A.  GURSON'S  VOID  MODEL 


Yielding  and  plastic  deformation  in  the  element  follows  the  model  proposed  by 
Gurson  for  symmetric  deformations  around  a  spherical  void  [17].  The  yielding  condition  is 


fxY 

K&0) 


2qj< Dcoshf-I^l  -  (l  +  q3<&2)=0 


(35) 


V  *  (JO  J 


where  <Pis  the  current  porosity,/*  is  the  hydrostatic  stress,  q  is  the  effective  stress,  and  cr0 
is  the  current  yield  stress.  This  model  assumes  equivalent  yield  stress  in  both  tension  and 
compression.  The  constants  qx,  q2,  and  q3  were  introduced  by  Tvergaard  in  order  to  provide 
a  better  match  with  numerical  studies  [16].  Aravas  provides  a  detailed  explanation  of 
implementing  this  model  in  a  static  finite  element  algorithm  for  three-dimensional  solid 
elements  [21].  The  procedure  used  here  is  essentially  the  same,  with  some  modifications  due 
to  the  different  element  formulation.  The  stress  is  then  transformed  to  local  coordinates,  and 
the  internal  force  vector  is  computed  as  described  above.  One  thing  to  be  noted  in  Eq.  (35) 
is  that  if  0  is  initially  0,  the  yielding  criteria  surface  is  identical  to  the  von  Mises  yield 
condition. 

After  calculating  the  strain  tensor  using  Eq.  (23),  any  previous  plastic  strain  is 
subtracted  using 

(36) 
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The  stress  is  then  calculated  using  Eq.  (24)  with  the  components  of  {se}-  Using  these 
values,  the  hydrostatic  stress,  deviatoric  stress  and  effective  stress  are  calculated  as  follows: 


P  =  -1A{<7XX 

(37) 

M=M  +  p{#} 

(38) 

y  =  ^i{s2l  +  S2  +  S3  +  2(s*  +  Ss  +  S6  )) 

(39) 

where  8  is  the  Kronecker  delta  function: 

{<?}  =  {lllOO}r  (40) 

At  this  point,  F  is  calculated  from  Eq.  (35).  If  F  is  greater  than  zero,  indicating  plastic  flow, 
iteration  is  required  to  determine  the  new  porosity  and  change  in  plastic  strain.  The 
predictor-corrector  method  used  in  Aravas  [21]  is  also  used  here.  Using  the  values  ofp  and 
q  previously  calculated  as  a  first  guess,  correction  factors  are  calculated  using 

{C,MA]-'{A1}  .(41) 

where 


-A 


dF 


£p  dq 


_  Ao 

^ Oq  Bn 


(42) 


w- 


f  +Asa 

aj  q 


K 


Kfd 


8F  d°o 


5p  dob  a tep 


Hj 

#F 

cpdob  dtep. 

+A& 


&F  foo  aF 


P  dfkj0  dtep  dp 


+As, 


•3G« 


\£Fj 

1  &F  foo 

dqda0  d/Ssq. 

dF  &r0 


aq  da0  dd£q 


'  dqdoo  dteq 


(43) 


and 


f =2q,® 


I $-2q,* 


-3<hP 


=  -2  qx<bp\ 


:-^r  +  2<?1<D^fsinh 


coshf  zlM.)  +  i^sinhl 


&L  =Zl 

dqdo-0  3 


The  values  for  a  and  p  are  used  to  correct  the  change  in  strain  caused  by  hydrostatic  pressure 
and  the  change  in  strain  caused  by  effective  stress  (See  Eq.  (46)),  which  are  then  used  to 
determine  the  change  in  the  plastic  strain  vector,  as  shown  in  Eq.  (47). 

Asn;w  =  As^d  +  a  (A 4=0) 

AsT=Asf  +  p  (A s°g  =  0)  (  } 


{a£'Ha*,{,5}+A*,U-){S} 

\2(l) 


This  change  in  plastic  strain,  { Asp  },  is  then  added  to  the  total  plastic  strain,  and  the  new 
elastic  strain  is  calculated  using  Eq.  (36),  and  the  stress  vector  is  calculated  again  using  Eq. 
(24).  Next,  the  change  in  void  content,  or  porosity,  is  calculated  as: 

A  growth  A  (Enucleation  (48) 


A4v-,=(i-®)2>;: 
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2 


(50) 


where  is  the  volume  fraction  of  void  nucleating  particles  and  sN  and  s*,  are  the  mean  and 
standard  deviation  of  a  normal  distribution  of  nucleation  strain,  as  suggested  by  Chu  and 
Needleman  [15]  and  utilized  by  Aravas  [21].  Then  the  effective  plastic  strain,  and  void 
content  are  updated  with 

<»> 

(7-0  )<jo 

sf+Al  =  sf  +  Asp  (52) 

where  sf  is  the  effective  plastic  strain  from  the  previous  time  step.  At  this  point,  the  change 
in  yield  stress  due  to  strain  hardening  is  calculated  (see  below).  If  either  a  or  /?  is  greater 
than  a  predetermined  tolerance,  the  process  iterates  beginning  with  Eq.  (41).  The  tolerance 
used  for  all  examples  presented  in  the  paper  is  1.0  x  10"4. 


B.  IMPROVING  STABILITY  IN  THE  CONSTITUTIVE  EQUATIONS 

Stability  in  the  procedure  outlined  above  is  very  dependent  on  the  order  in  which  the 
various  equations  are  evaluated.  Aravas  discusses  the  stability  problem  when  considering 
problems  involving  large  plastic  strains  [21].  While  the  change  in  plastic  strain  from  one 
time  step  to  the  next  will  theoretically  be  small,  high  strain  rate  and  changes  in  the  strain¬ 
hardening  characteristics  of  the  material  act  to  degrade  stability.  Equations  (42),  (43)  and 
(45)  are  not  complete  in  that  they  do  not  contain  all  of  the  derivative  terms  required  by  the 
chain  rule. 
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The  predictor-corrector  method  used  is  based  on  Newton's  method,  where  the 
correction  factors,  a  and  p,  are  found  from 


f\Mp  _  \f\ 

flM  f \p j  \fi 


where 


/,=  X  +29l(Dcosh  -^  -(l  +  ^O2) 


f2.As,SL  +  As& 

dq  dp 

f  ~Kd^  i  dA  50  i  da° 
1Mp  dp  dO  dAsp  dcr0  dAsp 

\  . 

f  -  3G  —  i  ^  —  i  ^  d(7° 
lMq  dq  dO  dAsq  da0  dAsq 


d2/i  dap  [  Aj_  Ktfj\  |  a2/,  do  [  a2/!  da, 

dcr0dq  dAs  q  dp2  dpdO  dAs  dpdcr0  dAs 


/,„  =|l+4^#-^^l  +  A£J-3G^+-^i^l  (59) 

’  q  dp  q]JdpdO  dAsq  dpdcr0  dAeqj  [  dq1  dqdcr0  dAsq 

Derivative  terms  that  are  zero  have  been  dropped,  and  K  and  G  are  the  bulk  and  shear 
moduli,  respectively,  f  is  Gurson's  void  function,  and^  is  the  flow  rule,  both  of  which  are 
driven  to  zero.  As  void  content  increases,  the  number  of  iterations  required  to  achieve 
convergence  increases  dramatically,  and  often  diverges,  instead.  By  removing  all  partial 
derivatives  relating  to  void  content,  the  stability  of  the  algorithm  is  greatly  increased,  at  the 


cost  of  only  2  to  3  extra  iterations,  depending  on  the  current  void  content.  The  dependence 
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on  the  current  yield  stress  is  retained  to  allow  convergence  across  a  transition  in  the  slope  of 
the  yield  stress  versus  strain  relationship. 

The  partial  derivatives  of  yield  stress  with  respect  to  ASp  and  Asq  are  computed  by 
calculating  the  slope  from  the  current  yield  stress  to  the  yield  stress  corresponding  to  the 
increase  in  plastic  strain  from  the  corrected  values  of  Asp  and  Asq.  This  means  that  these 
derivatives  will  be  zero  on  the  first  iteration,  since  both  of  the  control  variables  are  initialized 
to  zero.  Several  error  traps  must  be  included  to  prevent  round-off  and  truncation  error  from 
causing  the  algorithm  to  diverge,  and  to  prevent  porosity  from  decreasing  or  becoming 
negative.  One  of  the  assumptions  used  in  this  implementation  is  that  once  voids  form,  they 
do  not  disappear.  In  other  words,  voids  do  not  disappear  when  the  element  is  placed  in 
compression. 

C.  STRAIN  HARDENING 

Modeling  the  nonlinear  elasto-plastic  behavior  of  the  material  used  is  simplified  by 
constructing  a  piece-wise  linear  version  of  the  stress-strain  plot.  Using  the  tangent  modulus, 
Et,  for  each  piece- wise  region,  the  yield  stress  is  calculated  with 

+  'YjErAei  -  £1-1 )  +  Eu{ef  -  Sj-i )  (60) 

i~J 

where  a'o^  is  the  yield  stress  for  this  time  step,  cr°  is  the  original  yield  stress,  ef  is  the 
total  effective  strain  for  the  time  step  under  consideration,  and  s,  is  the  upper  strain  limit  of 
the  rth  linear  segment,  so  is  the  original  yield  strain,  and  is  calculated  by  the  program  as 
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simply  <Jq/E  .  The  current  effective  elastic  strain  is  calculated  by  dividing  the  current  yield 
stress  by  the  elastic  modulus,  E.  This  is  added  to  the  current  effective  plastic  strain  to  obtain 
the  total  effective  strain. 

As  an  example,  let  the  total  effective  strain  from  Eq.  (49)  lie  within  the  second  work¬ 
hardening  segment.  The  new  yield  stress  is 

cr‘o+A>  =  CTq  +  Eti  (  £1  -  So  )  +  En  (  £%td  ~  £1)  (61) 

Figure  3  illustrates  this  example.  The  algorithm  calculates  the  new  yield  stress  as  a  function 
of  the  cumulative  effective  plastic  strain,  the  current  effective  elastic  strain,  and  the  original 
yield  stress  for  each  iteration  of  the  damage  constitutive  equations. 


Figure  3.  Calculating  a  New  Yield  Stress. 
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IV.  HOURGLASS  MODE  CONTROL 


Only  one  integration  point  is  used  in  each  plane  parallel  to  the  reference  plane,  which 
results  in  under-integration  in  the  £-7  plane.  This  leads  to  spurious,  hourglass,  or  zero- 
energy  modes  in  the  element,  which  will  yield  useless  results  if  left  uncontrolled.  The  effects 
of  these  modes  are  shown  in  Fig.  4,  a  pinched  cylinder  where  one  eighth  of  the  structure  is 
modeled  by  utilizing  symmetry  boundary  conditions.  Clearly,  no  useful  information  can  be 
obtained  from  the  resulting  solution. 


Figure  4.  Hourglass  Modes  in  a  Pinched  Cylinder  Model. 
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Belytschko,  et  al.,  proposed  an  efficient  means  of  controlling  the  hourglass  modes 
of  a  similar  element  [25].  The  method  described  uses  a  portion  of  a  stiffness  matrix 
generated  by  full  integration  in  all  directions  to  modify  the  stiffness  matrix  generated  by 
under-integration.  Although  the  formulation  proposed  in  this  paper  does  not  use  a  stiffness 
matrix,  a  similar  approach  is  just  as  effective  in  controlling  these  modes. 

Rather  than  fully  integrating  in  all  three  directions,  the  element  is  fully  integrated  in 
the  £-77  plane,  but  under-integrated  in  the  ^ direction,  as  shown  in  Fig.  5.  The  procedure 
described  above  for  calculating  the  internal  force  vector  is  followed  to  generate  an  internal 
force  vector  related  to  these  four  integration  points. 


Figure  5.  Hourglass  Mode  Control  Integration  Points. 

The  algorithm  for  calculating  strain,  stress,  and  force  for  the  hourglass  mode  control 
integration  points  is  identical  to  the  algorithm  used  to  produce  the  main  internal  force  vector, 
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with  the  exception  of  plastic  strain.  The  damage  constitutive  equations  described  in  the 
previous  section  are  not  utilized  for  hourglass  mode  control.  The  internal  forces  generated 
from  the  two  integration  schemes  are  treated  like  the  stiffness  matrices  in  Belytschko  et  al. 
[25].  For  nz  integration  points  through  the  element  thickness,  the  new  force  vector  is 

{/j- {/£'■■’ }+*{/_}  («> 

where 

{/_}-{/£“ }-{/“"}  («3) 

and 

f  t 2 

h=—  (64) 

A 

The  variable  h  is  used  here  in  Eqs.  (62)  and  (64)  instead  of  the  fused  in  Belytschko 
et  al.  [25],  to  avoid  confusion  with  the  various  strains  discussed  in  this  paper.  The  variables 
used  to  calculate  h  are  the  element  thickness  (t)  and  the  element’s  surface  area  (A).  The 
effect  of  r  follows  that  described  in  Belytschko  et  al.  [25],  and  is  set  to  0.05.  The  range  of 
values  for  r  that  effectively  controls  the  hourglass  modes,  but  does  not  greatly  affect  the 
overall  element  stiffness  is  roughly  0.046  to  0.057  (determined  experimentally).  Since  the 
elements  are  in  arbitrary  orientation  in  3-D  space,  the  area  calculation  is  computed  as  the 
sum  of  the  area  of  the  two  triangles  formed  by  dividing  the  element  at  the  diagonal  between 
nodes  one  and  three: 

-  D/  +  lU24  -  D\  (65) 

where 
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Di  =  (x\  -  xf  )  (x\ -x\)  +  (x\  -  x\  )  (x\ -x\)  +  (x\  -  x]  )  (xl  -  X* ) 


D2  =  (x\  -  xl )  (x\  -xl)  +  (x\  -  xl )  (x\  -xl)  +  (x\  -  xl )  (xl  -  xl ) 


(66) 


and 


h  =  j(xl-xl)2  +  (x22  -x\f  +  (xl  -x])2 
l2  =  ^(xl-xl)2  +  (xl-xl)2  +  (xl-xl)2 


(67) 

W(*i  -xf)2  +  (xl-xl)2  +  (xl  -xl)2 

u=Ax\  -*i4)2  +  ix\  ~xl)2  +  (xl  -  4  )2 

and  (xf  ,xl,xl)  is  the  location  of  node  k  in  the  global  coordinate  system.  For  these 

calculations,  the  element  is  assumed  to  be  flat  (no  curvature  along  either  the  £  or  rj 
directions).  The  effectiveness  of  this  method  of  control  is  shown  in  Fig.  6,  the  same  problem 
as  shown  in  Fig.  4,  but  with  the  hourglass  mode  control  described  above. 


Figure  6.  Pinched  Cylinder  Model  with  Hourglass  Mode  Control. 
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All  verification  problems  analyzed  in  the  following  section  were  completed  with 
hourglass  control  enabled.  The  current  implementation  uses  hourglass  control  consistently, 
but  allows  easy  modification  to  make  hourglass  control  a  user-defined  option.  The  internal 
hourglass  mode  control  can  be  disabled  when  the  element  is  used  in  a  general  finite  element 
program  that  includes  various  methods  of  spurious  mode  control. 
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V.  IMPLEMENTATION  IN  RESEARCH-ORIENTED  FE  CODE 


An  in-house  general  purpose  finite  element  analysis  program,  referred  to  herein  as 
FEA,  was  used  to  verify  the  algorithms  discussed.  Several  of  these  verification  problems  are 
presented  in  the  next  section.  The  program  is  written  in  FORTRAN,  and  is  tailored  for  easy 
modification,  rather  than  speed  or  ability  to  handle  very  large  problems.  The  program  reads 
the  system  parameters  from  a  single  input  file,  then  writes  the  results  to  two  ASCII  text  files: 
one  contains  the  nodal  coordinates  and  displacements  for  all  time  steps,  the  other  contains 
the  data  associated  with  each  element  for  all  time  steps  (stress  and  strain  tensors,  effective 
plastic  strain,  void  content,  yield  stress,  etc.). 

A.  PREPROCESSOR 

Since  this  program  is  locally  generated,  there  is  no  preprocessor  or  postprocessor 
available.  In  order  to  facilitate  running  several  problems,  a  preprocessor  was  written  in 
MATLAB  that  generates  the  proper  input  file  for  several  simple  geometric  shapes.  The 
MATLAB  script  file  will  create  a  properly  meshed  flat  plate,  open  box,  cylinder,  complete 
spherical  cap,  or  a  spherical  cap  with  a  hole  in  its  top.  The  listing  for  the  preprocessor  is 
presented  in  Appendix  A. 

B.  POSTPROCESSOR 

A  postprocessor,  also  written  as  a  MATLAB  script  file,  utilizes  a  graphical  user 
interface  (GUI)  to  create  the  required  graphs.  It  also  has  the  capability  to  create  an  animated 


31 


display  of  the  displacement  over  time.  The  displacements  are  scaled  to  make  them  clearly 
visible.  This  feature  is  useful  in  determining  the  effectiveness  of  hourglass  control  and 
verifying  the  boundary  conditions  applied.  However,  the  MATLAB  movie  is  memory 
intensive  and  slow,  and  is  therefore  useful  only  as  a  qualitative  verification  of  proper 
problem  definition.  All  graphs  of  problem  solutions,  except  as  noted,  were  generated  using 
this  postprocessor.  The  full  listing  of  all  associated  files  is  provided  in  Appendix  B. 

C.  FEA  SHELL  ELEMENT  SUBROUTINE  DESCRIPTION 

The  FORTRAN  source  code  for  the  FEA  implementation  of  the  element  presented 
here,  along  with  the  proposed  hourglass  control  and  constitutive  model,  is  listed  without 
comment  in  Appendix  C.  Since  the  primary  purpose  of  the  FEA  program  is  to  assist  in  the 
development  of  finite  element  algorithms,  this  subroutine  is  written  with  readability,  rather 
than  computational  efficiency,  as  a  primary  concern,.  A  brief  outline  of  the  subroutine  is 
shown  in  Fig.  7.  The  FEA  program  treats  each  type  of  element  (shell,  beam,  etc.)  in  a  single 
block,  and  loops  through  each  element,  then  each  integration  point  within  the  element. 

The  order  of  calculation  when  evaluating  the  constitutive  equations  is  critical,  and 
any  variation  will  cause  instability.  The  two  control  variables,  dep  and  deq,  are  updated  by 
the  correction  factors  on  each  iteration.  All  other  variables  are  initialized  to  the  values  from 
the  previous  time  step  within  each  iteration.  This  essential  feature  allows  stability  of  the 
algorithm.  As  discussed  above,  the  variation  of  void  content  with  dep  and  deq  is  not  used. 
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but  the  variation  of  yield  stress  with  the  two  control  variables  is  calculated  and  included  in 


the  full  expansion  of  the  partial  derivatives. 
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Program  Outline  for  ELSH4L.F 

1 .  Enter  subroutine 

2.  Set  constants 

3 .  Load  gauss  points  and  weights 

4.  Begin  loop  over  each  shell  element 

A.  Initialize  element  internal  force  vectors 

B.  Retrieve  element  material  properties 

C.  Compute  element  indices  into  system  vectors 

D.  Load  nodal  coordinates  and  displacements 

E.  Compute  local  coordinate  system  and  surface  area 

F.  Compute  rotation  matrices 

G.  Transform  displacements  to  local  coordinate  system 

H.  Calculate  hourglass  force  vector 

I.  Begin  Integration  Loop 

1)  Load  shape  functions  and  derivatives 

2)  Calculate  Jacobian,  its  inverse,  and  its  determinate 

3)  Calculate  strain-displacement  matrix  ([B]) 

4)  Compute  strain  ({s}  =  [B]  {d}) 

5)  Transform  strain  to  local  coordinate  system 

6)  Subtract  any  previous  plastic  strain  from  strain  tensor 

7)  Load  last  values  of  o0,  <|>,  and  spdT 

8)  Compute  stress  ({a}  =  [D]  {s}) 

9)  Compute  yield  function  (F) 

10)  If  F  >  0,  set  dep  =  deq  =  0  and  begin  iteration 

a)  Calculate  required  partial  derivatives 

b)  Prevent  Numerical  Errors  (Divide  by  Zero,  etc.) 

c)  Retrieve  Previous  <|>,  eV ,  and  cj0 

d)  Update  dep  and  deq 

e)  Compute  Incremental  Plastic  Strain  Tensor 

f)  Compute  Void  Growth 

g)  Compute  New  cry,  Sy,  and  spefr 

h)  Compute  Void  Nucleation  Factor 

i)  Compute  New  <|)  and  aQ 

j)  Recalculate  Yield  Function  (F) 

h)  If  F  is  not  Zero,  go  to  step  a) 

1 1)  Last  values  used  are  retained. 

12)  Transform  stress  to  global  coordinate  system 

13)  Calculate  internal  force  for  this  integration  point  and  sum 

1 4)  Store  material  tensors  and  constants  in  system  vectors 

15)  End  of  integration  loop 

J.  Apply  hourglass  correction  force 

K.  Transform  internal  force  vector  to  global  coordinate  system 

L.  Compute  element  mass  vector 

M.  Load  mass  and  force  vectors  into  system  vector 

N.  End  of  element  loop 

5.  End  of  subroutine 


Figure  7.  Outline  of  FEA  Implementation. 
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VL  NUMERICAL  EXAMPLES 


The  element  presented  here  has  been  extensively  tested  using  a  variety  of  problems 
in  which  an  analytic  solution  was  available,  and  has  produced  satisfactory  results  in  all  cases 
studied.  The  transformation  matrices  and  constitutive  equations  were  verified  with 
specifically  tailored  problems,  and  the  element  passes  the  patch  test.  The  examples  presented 
below  are  representative  of  the  test  cases  used  to  verify  the  element,  and  are  presented  in  the 
order  of  curvature:  flat,  singly  curved,  and  doubly  curved.  For  each  case,  an  elastic  solution 
is  compared  to  available  results,  and  then  the  elasto-plastic  solutions  are  presented. 

A.  ELASTIC  PLATE 

A  plate  clamped  on  all  four  sides  is  subjected  to  a  concentrated  force  at  its  center. 
The  elastic  modulus  is  10  Msi  (68.95  GPa),  the  density  is  0.1647  slugs/in3  (0.147  kg/cm3), 
and  Poisson’s  ratio  is  0.2.  The  dimensions  of  the  plate  are  10  in  x  10  in  x  0.1  in  thick 
(25.4cm  x  25.4  cm  x  0.254  cm).  The  yield  stress  is  set  high  enough  to  ensure  a  completely 
elastic  response.  Two  integration  points  through  the  thickness  are  used.  The  applied  force 
is  40  lbf  (177.9  N).  The  finite  element  mesh  uses  symmetry  to  model  one  quarter  of  the 
plate,  with  appropriate  symmetry  boundary  conditions  applied.  Both  a  2x2  (4  element)  and 
4x4  (16  element)  mesh  are  used  in  the  finite  element  analysis.  The  results  for  both  meshes 
and  the  analytic  solution  are  shown  in  Table  1 .  Using  a  four-element  mesh,  the  new  shell 
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element  obtained  a  displacement  within  3.27%  of  the  analytic  solution,  and  0.82%  using  a 
16  element  mesh. 


Table  1 .  Comparison  of  Results  for  Elastic  Clamped  Plate. 


Analysis  Type 

Center  Node  Peak  Displacement  (in) 

2  x  2  FE  Mesh  (Dynamic) 

-4.74xl0'2 

4  x  4  FE  Mesh  (Dynamic) 

-4.94x1  O'2 

Analytic  (Twice  the  Static  Solution) 

-4.90xl0-2 

B.  THICK  CLAMPED  PLATE  UNDER  PRESSURE  LOAD  IN  ELASTO- 
PLASTIC  REGION 

A  thick  steel  plate,  clamped  on  all  four  sides,  is  subjected  to  dynamic  pressure  load. 
The  plate  is  6m  by  6m  by  0.6m  thick  (for  a  10:1  ratio).  Table  2  shows  the  material 
properties  for  the  structure.  Table  3  shows  the  properties  for  the  void  model. 


Table  2.  Material  Properties  of  Clamped  Plate. 


Property 

Value 

Units 

Elastic  Modulus  (E) 

2xl0u 

Pa 

Tangent  Modulus  (Ej.) 

2xi010 

Pa 

Density  (p) 

7850 

kg/m3 

Poisson’s  ratio  (v) 

0.29 

(none) 

Yield  Stress  (Syp) 

2.5x10s 

Pa 
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Table  3.  Void  Characteristics  of  Clamped  Plate. 


Initial  Void  Content  (<t>0) 

0.0 

Nucleating  Particle  Content  (<DN) 

0.04 

Mean  Nucleation  Strain  (e^ 

0.3 

Nucleation  Strain  Standard  Deviation  (%) 

0.1 

Model  Constant  q, 

1.5 

Model  Constant  c ^ 

1.0 

Model  Constant  ^  (  =  qi2) 

2.25 

One  quarter  of  the  plate  is  modeled  with  a  three  by  three  element  mesh,  for  a  total  of 
nine  elements,  with  appropriate  symmetry  boundary  conditions  applied,  and  is  shown  below 
in  Fig.  8. 


Figure  8.  Nine  Element  Clamped  Plate  Mesh. 
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Four  integration  points  are  used  through  the  thickness  of  each  element.  The 
calculation  time  step  is  10'5  seconds,  while  the  data  is  plotted  at  2x1  O'4  second  increments. 
The  analysis  terminates  at  0.04  seconds.  The  plate  is  subjected  to  a  uniformly  distributed 
pressure  acting  downwards,  and  includes  the  appropriate  drilling  moment.  The  pressure 
increases  linearly  from  0.0  Pa  at  the  start  to  80  MPa  at  0.01  seconds,  then  remains  constant 
for  the  duration  of  the  analysis.  Three  cases  were  analyzed:  no  void  effects,  void  growth 
effects  only,  and  void  growth  and  nucleation.  Figure  9  illustrates  the  effect  of  voids  on  the 
effective  stress  versus  the  effective  strain  relationship  in  the  top  of  one  of  the  border 
elements.  When  void  effects  are  not  included,  the  Von-Mises  Equivalent  (VME)  stress 
follows  the  yield  stress  in  the  plastic  region  (See  Fig.  9a). 

Including  void  effects  causes  the  yielding  before  the  VME  stress  reaches  the  yield 
stress.  As  the  void  content  increases  with  continued  plastic  flow,  the  VME  stress  falls  farther 
below  the  yield  stress,  as  shown  in  Fig.  9b  (see  Eq.  (35)).  Table  4  summarizes  the  results 
of  the  three  cases.  The  most  significant  difference  is  in  the  transverse  normal  stress,  cte. 
This  difference  resulted  in  a  2.0%  increase  in  the  effective  plastic  strain,  and  a  2.6%  increase 
in  the  peak  deflection  of  the  center  of  the  plate. 

Another  interesting  point  is  that  the  transverse  normal  stress  was  compressive  and 
constant  throughout  the  thickness  up  until  plastic  flow,  as  assumed  in  the  formulation,  then 
varied  as  the  stress  in  the  bottom  fiber  decreased  and  became  tensile.  Figure  10  shows  the 
variation  of  transverse  normal  stress  through  the  thickness  of  the  element  at  the  time  of  peak 
stress.  Eventually,  the  transverse  normal  stress  varied  from  compressive  at  the  top,  where 
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the  pressure  was  applied,  and  tensile  at  the  bottom.  The  examples  that  follow  will  not 
include  a  separate  analysis  for  void  growth  effects  only,  since  the  difference  of  including 
nucleation  effects  at  the  resulting  small  plastic  strains  obtained  do  not  cause  significantly 
different  results. 
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Figure  9.  Void  Effects  in  a  Clamped  Plate  with  Pressure  Loading:  Top  (a)  -  No  Void 
Effects,  Bottom  (b)  -  Void  Effects. 
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Table  4.  Summary  of  Results  for  Clamped  Plate  with  Pressure  Load. 


Peak  Values  for  Element  #3 

No  Void  Effects 

Void  Growth 

Growth  and  Nucleation 

<?VM  (GPa) 

1.4498 

1.3789 

1.3784 

^effective 

0.0588 

0.0597 

0.0597 

e 

1.4501 

1.3603 

1.3598 

a™  (GPa) 

0.6561 

0.6030 

0.6073 

(GPa) 

(Max  Tension) 

-0.0003 

0.0065 

0.0091 

<7=  (GPa) 

(Max  Compression) 

-0.0363 

-0.0322 

-0.0320 

^plastic  (effective) 

0.0540 

0.0551 

0.0551 

4>  (Porosity) 

NA 

0.0355 

0.0357 

Deflection  (m) 

(Center  Node) 

-0.3801 

-0.3897 

-0.3898 

C.  THICK  CLAMPED  PLATE  WITH  CENTRAL  POINT  LOAD  IN  THE 

ELASTO-PLASTIC  REGION 

The  geometry  and  material  properties  of  this  example  are  identical  to  those  used  in 
the  previous  example.  The  pressure  load  has  been  replaced  with  a  point  load  at  the  center 
node.  Four  cases  are  studied:  without  void  effects  or  a  drilling  moment,  with  void  effects, 
with  a  drilling  moment  applied,  and  with  both  void  effects  and  a  drilling  moment  applied. 
All  four  sides  are  clamped,  and  the  center  node  displacement  is  restricted  in  the  x  and  y 
directions.  Then,  a  DM  equal  to  the  applied  force  times  one-half  the  plate  thickness  is 
applied  for  both  the  no-void  and  void  cases.  The  time  history  of  the  stress  components  in  the 
center  element  is  shown  in  Fig.  11,  and  the  relationship  between  porosity  and  effective 
plastic  strain  is  illustrated  in  Fig.  12.  The  results  for  all  six  cases  are  summarized  below  in 
Table  5. 
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Figure  1 1 .  Stress  Component  Time  History  in  Bottom  Fiber:  Clamped  Plate  with 
Concentrated  Load,  Void  Effects,  and  Drilling  Moment  Applied. 
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Table  5.  Summary  of  Results  for  Clamped  Plate  Subjected  to  Point  Force. 


Peak  Values  for  Center 
Element 

No  Voids 

No  DM 

Voids 

No  DM 

No  Voids 
Drilling  Moment 

Voids 

Drilling  Moment 

°vm  (GPa) 

2.1119 

1.9021 

2.0540 

1.8590 

^effective 

0.0892 

0.0929 

0.0865 

0.0887 

a„  (GPa) 

2.0008 

1.8407 

2.0000 

1.7498 

aw  (GPa) 

1.9588 

1.7245 

1.8589 

1.7234 

an  (GPa) 

(Max  Tension) 

0.0048 

0.0213 

-0.0003 

-0.0003 

azz  (GPa) 

(Max  Compression) 

-0.0036 

-0.0057 

-0.1751 

-0.1691 

gptoic  (effective) 

0.0838 

0.0880 

0.0812 

0.0840 

4>  (Porosity) 

NA 

0.0663 

NA 

0.0601 

Deflection  (m) 

(Center  Node) 

-0.4413 

-0.4544 

-0.4350 

-0.4451 

The  drilling  moment  increases  the  effective  stress  on  the  fiber  in  compression,  and 
decreases  the  effective  stress  on  the  fiber  in  tension.  The  reduction  in  the  tensile  stress 
results  in  a  reduction  in  the  effective  plastic  strain.  Including  void  effects  decreases  the 
VME  stress  on  the  fiber  in  tension,  and  slightly  increases  the  VME  stress  on  the  fiber  in 
compression.  The  effective  plastic  strain  is  also  increased.  All  of  these  results  indicate  that 
this  formulation  correctly  predicts,  in  a  qualitative  sense,  the  effects  of  drilling  moments  and 
void  growth  and  nucleation.  Including  only  void  effects  increased  the  effective  plastic  strain 
by  5.0%,  while  including  only  the  drilling  moment  decreased  the  effective  plastic  strain  by 
3.1%.  Including  both  void  effects  and  the  drilling  moment  increased  the  effective  plastic 
strain  by  0.2%,  indicating  the  importance  of  applying  both  effects  together. 
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D.  SIMPLY  SUPPORTED  PLATE  WITH  CENTRAL  POINT  LOAD  IN  THE 
ELASTO-PLASTIC  REGION 

The  material  properties,  geometry,  mesh,  and  time  values  from  the  clamped  plate 
problem  above  are  used  here.  The  magnitude  of  the  force  is  8x10*  N,  and  the  drilling 
moment,  when  applied,  is  -2.4x10s  N.  The  same  four  cases  are  analyzed:  no  void  or 
drilling  moment  effects,  drilling  moment  effects  only,  void  effects  only,  and  both  void 
and  drilling  moment  effects.  The  analysis  results  for  all  four  cases  are  summarized  in 
Table  6. 


Table  6.  Summary  of  Results  for  a  Simply  Supported  Plate  with  a  Point  Force. 


Peak  Values  for  Center 
Element 

No  Voids 
No  DM 

No  Voids 
Drilling  Moment 

Voids 

No  DM 

Voids  and 
Drilling  Moment 

'c? 

1 

! 

3.0874 

3.0089 

2.5061 

2.4881 

0.1360 

0.1321 

0.1450 

0.1393 

(GPa) 

3.0253 

2.9123 

2.4686 

2.4206 

(GPa) 

3.0459 

2.8889 

2.4612 

2.3928 

(GPa) 

(Max  Tension) 

0.0060 

-0.0003 

0.0114 

-0.0003 

tfzz  (GPa) 

(Max  Compression) 

-0.0035 

-0.1802 

-0.0035 

-0.1722 

gpiastic  (effective) 

0.1242 

0.1382 

0.1327 

(Porosity) 

NA 

NA 

0.1197 

0.1092 

Deflection  (m) 

(Center  Node) 

-0.9219 

-0.8990 

-0.9770 

-0.9420 

The  qualitative  results  for  this  example  are  die  same  as  for  the  clamped  plate.  Since 
the  applied  force  causes  greater  initial  yielding,  void  growth  and  nucleation  has  a  greater 
effect.  The  drilling  moment,  when  added  to  the  void  effects,  reduces  the  amount  of  effective 
plastic  strain  in  tension  and  displacement. 
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E.  ELASTIC  PINCHED  CYLINDER 


An  open-ended  cylinder  of  radius  5.0  in.  (12.7  cm),  length  10.35  in.  (26.289  cm),  and 
thickness  0.094  in.  (0.23876  cm)  is  subjected  to  a  pinching  load  of  100  lbf  (444.82  N)  (See 
Fig.  6).  The  elastic  modulus  is  10.5  msi  (72.4  Gpa),  Poisson’s  ratio  is  0.3125,  and  the 
density  is  3.125xl0'3  slugs/in3  (2783  kg/m3).  The  load  is  applied  as  a  step  function  beginning 
at  time  t  =  0  seconds.  Using  symmetry,  the  problem  was  reduced  to  a  one-eighth  section  of 
the  cylinder.  The  dynamic  value  should  be  twice  the  analytic  static  value.  Inextensional 
shell  theory  gives  a  static  radial  contraction  of  0.1117  in.  (0.28372  cm).  The  maximum 
radial  contraction  of  the  model  is  0.1995  in.  (0.5067  cm),  which  translates  to  a  static  radial 
contraction  of 0.09975  in  (0.2534  cm).  Using  a  256-element  mesh  (16  by  16),  the  maximum 
radial  contraction  was  0.2207  in.  (0.5606  cm),  for  a  static  contraction  of  0.1 104  in  (0.2804 
cm).  The  results  are  summarized  in  Table  7.  The  16-element  solution  is  within  10.7%  of 
the  analytic  solution,  while  the  256-element  mesh  is  within  1.21%. 


Table  7.  Comparison  of  Results  for  Elastic  Pinched  Cylinder. 


Analysis  Type 

Radial  Contraction  (in) 

Finite  Element  with  16  Element  Mesh 

0.1995  (0.5067  cm) 

Finite  Element  with  256  Element  Mesh 

0.2207  (0.5606  cm) 

Analytic  (Twice  the  Static  Solution) 

0.2234  (0.5674  cm) 

F.  THICK  PINCHED  CYLINDER  IN  THE  ELASTO-PLASTIC  REGION 

The  material  and  void  properties  shown  in  Tables  3  and  4  are  used  for  an  open-ended 
cylinder  with  a  pinching  force  at  its  center.  The  top  half  of  the  cylinder  is  modeled  using  a 
36  element  mesh  with  appropriate  symmetry  boundary  conditions  along  the  bottom  edge  of 
mesh,  as  shown  in  Fig.  13.  The  analysis  is  calculated  in  10'5  second  steps,  with  output  every 
2x1 04  seconds,  from  0.0  seconds  to  0.04  seconds.  The  pinching  force  is  66.792  kN  on  each 
side,  with  a  drilling  moment  of  848.258  Nm  applied  on  the  cases  indicated.  The  same  four 
cases  are  compared:  no  voids  or  drilling  moment,  void  effects  only,  drilling  moment  only, 
and  both  void  effects  and  drilling  moment  (DM).  The  effects  of  void  growth  and  nucleation 
and  the  drilling  moment  is  shown  in  the  stress  component  time  histories;  Fig.  14  shows  the 
stresses  without  voids  or  a  drilling  moment,  and  Fig.  15  shows  the  stresses  with  both  effects 
included.  The  results  for  all  cases  are  summarized  in  Table  8. 
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Table  8.  Summary  of  Results  for  Elasto-Plastic  Pinched  Cylinder. 


Peak  Values  for  Center 
Element 

No  Voids 

No  DM 

No  Voids 
Drilling  Moment 

Voids 

No  DM 

Voids  and 
Drilling  Moment 

avM  (MPa) 

262.40 

265.05 

259.28 

283.49 

1.6504 

1.5722 

1.4273 

2.7049 

(MPa) 

182.46 

160.38 

170.12 

150.31 

Gyy  (MPa) 

291.71 

283.77 

291.83 

297.97 

csa  (MPa) 

(Max  Tension) 

5.0619 

-0.0238 

1.0270 

7.5235 

oa  (MPa) 

(Max  Compression) 

-28.739 

-7.0364 

-35.373 

gpiasic  (effective)  (x  103) 

0.7828 

0.7583 

0.5771 

1.6923 

4>  (Porosity)  (x  103) 

NA 

NA 

0.4773 

1.1970 

Deflection  (mm) 

(Center  Node) 

Global  Maximum 

2.7589 

12.919 

-0.0220 

-0.0214 

Deflection  (mm) 

(Center  Node) 

Global  Minimum 

-31.750 

-31.750 

-31.750 

-31.750 

Since  there  is  only  a  small  amount  of  plastic  strain,  the  effects  of  void  growth  and 
nucleation  are  greatly  reduced.  Even  with  this  small  amount  of  plastic  flow,  die  effective 
plastic  strain  was  increased  by  over  116  percent  with  the  inclusion  of  both  drilling  moment 
and  void  effects,  while  the  peak  center  node  displacement  was  only  increased  by 
approximately  4.5  percent. 
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Figure  13.  Mesh  Structure  for  Pinched  Cylinder. 


Figure  14.  Stress  Component  Time  History  in  Bottom  Fiber:  Pinched  Cylinder. 
No  Void  Effects  or  Drilling  Moment. 


Figure  15.  Stress  Component  Time  History  in  Bottom  Fiber:  Pinched  Cylinder,  Void 
Effects  and  Drilling  Moment  Applied. 


G.  ELASTIC  SPHERICAL  CAP  WITH  A  CENTER  HOLE 


The  results  of  analysis  with  the  shell  element  developed  here  are  compared  to  results 
of  the  problem  proposed  in  MacNeal  and  Harder  [26].  The  structure  is  a  hemisphere  of 
radius  1 0  units  with  a  thickness  of  0.04  units,  and  has  an  1 8°  hole  cut  in  the  center.  Taking 
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advantage  of  the  symmetry  in  the  structure,  only  %  of  the  structure  is  modeled.  An  eight-by¬ 
eight  mesh  is  used,  for  a  total  of  64  elements,  with  four  integration  points  through  the 
thickness  of  each  element.  The  mesh  used  is  shown  in  Fig.  16.  The  material  properties  of 
the  structure  are:  E  =  6.825xl07,  p  =  0.001,  and  v  =  0.3.  No  void  or  drilling  moment  effects 
are  employed  in  this  example.  Opposing  forces  of  magnitude  2.0  are  applied  at  each 
quadrant:  the  force  at  node  7,3  is  of  magnitude  -1 .0  and  parallel  to  the  y-axis,  and  the  force 
at  node  one  is  of  magnitude  1.0  and  parallel  to  the  x-axis.  To  obtain  a  representative  static 
response  a  dynamic  model,  the  applied  force  begins  with  0.0  magnitude  at  time  0,  and 
increases  linearly  with  a  rise  time  of  0. 1 6  seconds  to  the  specified  value.  A  calculation  time 
step  of  lxl O'6  seconds  is  used,  with  termination  at  0.22  seconds.  The  maximum  deflection 
of  node  one  was  0.0929,  where  the  theoretical  value  is  0.0940,  for  a  normalized  displacement 
of  0.988.  This  is  within  the  range  of  the  results  listed  in  MacNeal  and  Harder  from  the 
QUAD2  and  QUAD4  elements,  which  are  considered  to  be  accurate  for  this  problem  [26]. 
These  results  are  summarized  in  Table  9. 


Table  9.  Comparison  of  Results  for  Spherical  Cap  with  Elastic  Loading,  and  Results  from 
MacNeal  and  Harder  [26]. 


New  Shell  Element 

QUAD2 

QUAD4 

Normalized  Displacement 

0.988 

0.986 

1.005 
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Figure  16.  Mesh  Structure  for  Spherical  Cap. 

H.  SPHERICAL  CAP  WITH  A  CENTER  HOLE,  IMPACT  LOADING 

Key  and  Hoff  [1 1]  used  the  previous  problem  to  test  their  element  with  an  impact 
(step)  load.  The  structure  and  matrial  properties  are  the  same,  and  the  mesh  is  the  same  as 
shown  in  Fig.  1 6.  The  load  is  applied  at  its  maximum  at  t=0.0+  seconds,  rather  than  ramped 
up.  Figure  17,  which  plots  the  displacement  of  node  one  obtained  from  both  the  shell 
element  presented  here  and  the  element  developed  by  Key  and  Hoff  [11],  shows  that  the 
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Mean  Node  Displacement 


results  of  this  element  are  comparable  to  those  obtained  by  other  elements  under  impact 
loading. 


Figure  17.  Node  One  Displacement  Time  History:  Pinched  Spherical  Cap  with  Impact 
Loading. 

I.  SPHERICAL  CAP  WITH  A  CENTER  HOLE,  ELASTO-PLASTIC 
LOADING 

The  spherical  cap  structure  shown  in  Fig.  16  is  now  used  to  verify  the  stability  and 
convergence  of  the  damage-constitutive  equations  under  double  curvature,  and  to  illustrate 
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both  the  effects  of  void  growth  and  nucleation  and  a  drilling  moment  on  a  more  complex 
structure.  The  material  properties  used  are  the  same  as  shown  in  Tables  3  and  4.  The 
structure  has  a  radius  of  5.0  meters  and  a  thickness  of  25  cm,  for  a  radius  to  thickness  ratio 
of  20: 1 ,  commonly  considered  the  limit  of  thin-shell  theory.  Although  a  force  is  applied  at 
each  quadrant,  all  four  loads  are  directed  inwards.  Therefore,  only  one  load  is  applied  at 
node  37  of  the  mesh  shown  in  Fig.  16.  The  applied  load  is  5.9397xl06  N,  with  a  drilling 
moment  of  -7.4246x1 05  Nm  applied  for  the  cases  indicated.  A  calculation  time  step  of 
lxl  O'5  seconds  is  used,  with  output  every  5x10"*  seconds  and  stopping  at  0.2  seconds.  The 
load  is  initial  zero,  and  increases  linearly  to  its  maximum  at  0.01  seconds.  The  results  shown 
in  Table  10  are  for  the  element  nearest  the  applied  load,  where  stress  is  maximum.  The 
porosity  versus  effective  plastic  strain  relationship  is  shown  in  Fig.  18.  Due  to  the  small 
amount  of  plastic  strain,  porosity  is  restricted  to  the  linear  zone,  as  shown.  Figures  19a  and 
19b  illustrate  the  effective  stress  versus  effective  strain  in  the  same  element  on  the 
compressive  side  and  tensile  side,  respectively.  Note  that  as  the  structure  returns  from  a  peak 
displacement,  the  stress  follows  the  elastic  modulus,  not  the  tangent  modulus.  This  reflects 
the  correct  behavior  of  material:  the  elastic  modulus  is  not  affected  by  plastic  deformation, 
only  the  yield  stress  is  affected.  The  node  displacement  where  the  force  is  applied  is  shown 
in  Fig.  20.  This  figure  shows  that  there  are  at  least  two  vibration  modes  in  operation. 
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Table  10.  Summary  of  Results  for  Elasto-Plastic  Spherical  Cap. 


Peak  Values  for 
Element  #25 

No  Voids 

No  DM 

No  Voids 

Drilling  Moment 

Voids 

No  DM 

Voids  and 
Drilling  Moment 

Gvm  (MPa) 

250.26 

250.76 

25026 

250.75 

1.0730 

1.1021 

1.0731 

1.1021 

(MPa) 

202.09 

199.44 

202.08 

199.43 

CTw  (MPa) 

-47.093 

-58.183 

-47.087 

-58.192 

(j„  (MPa) 

-5.2570 

-12.433 

-5.2571 

-12.433 

2.1002 

4.1795 

2.1020 

4.1849 

<I>  (Porosity)  (x  106) 

NA 

NA 

6.7149 

12.763 

Deflection  (m) 

(Point  of  Loading) 

0.0233 

0.0246 

0.0233 

0.0246 

Effective  Stress  Effective  Stress 


Effective  Strain 


Effective  Strain 


Figure  19.  Comparison  of  Void  and  Drilling  Moment  Effects  in  a) 
Compression  (Top)  and  b)  Tension  (Bottom)  for  Spherical  Cap. 


57 


Node  #37 


Figure  20.  Contact  Node  Mean  Displacement  Time  History:  Spherical  Cap  with  Void  and  Drilling 
Moment  Effects. 

During  testing,  it  became  apparent  that  this  problem  is  not  suitable  for  testing  elasto- 
plastic  analysis.  Since  the  structure  is  concave,  the  entire  structure  quickly  collapses  once 
yielding  begins.  In  addition,  the  single  point  used  to  prevent  rigid-body  motion  also 
provided  a  stress  concentration  and  additional  yielding  (the  "comer"  began  folding  over). 

This  necessitated  using  a  force  that  just  causes  plastic  flow,  but  does  not  collapse  the 
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structure  or  cause  yielding  near  the  anchored  node.  This  is  illustrated  by  an  analysis  of  a  5 
percent  greater  load  than  used  for  the  analysis  shown  in  Fig.  20.  The  resulting  node  37 
displacement  is  shown  in  Fig.  21,  and  the  deformed  structure  in  Fig.  22.  However,  it  is  still 
apparent  that  the  qualitative  results  obtained  in  the  previous  elasto-plastic  examples  carried 
through  to  this  problem;  both  voids  and  the  drilling  moment  decreased  stress  in  fibers  under 
tension,  and  increased  stress  in  fibers  under  compression.  Due  to  the  small  amount  of  plastic 
strain,  the  increase  in  porosity  was  extremely  small,  which  minimized  the  effects  of  the 
voids. 


Time  (sec) 

Figure  21 .  Contact  Node  Displacement:  Spherical  Cap  with  5%  Greater  Load. 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 

A  new  shell  formulation  was  developed  for  transient  dynamic  analysis  which 
includes  both  void  effects  with  plastic  deformation  and  the  transverse  normal  stress  with 
drilling  moment.  The  element  is  compatible  with  most  shell  elements  which  have  three 
translation  and  three  rotation  degrees  of  freedom  per  node. 

A  numerically  stable  scheme  was  developed  for  Gurson's  nonlinear  constitutive 
equation  model.  The  model  includes  both  void  nucleation  and  void  growth.  Furthermore, 
hourglass  control  was  implemented  into  the  algorithm  to  avoid  spurious  modes  caused  by 
the  under-integration  scheme. 

The  drilling  moment  was  important  for  thick  plates  and  shells.  It  induced  large 
transverse  normal  stress  and  affected  the  plastic  deformation.  Similarly,  the  effect  of  voids 
on  plastic  deformation  was  not  negligible.  Numerical  studies  show  that  the  effective  plastic 
strain  increased  up  to  fifteen  percent  due  to  the  combined  effects  of  voids  and  the  transverse 
normal  stress. 

The  shell  element  presented  here  was  tested  on  many  cases,  some  of  which  had 
reference  solutions  for  comparison.  The  new  element  gave  accurate  results  to  those 
problems. 

Although  the  element  formulation  presented  here  peformed  well  in  numerical 
studies,  the  problems  studied  are  not  in  the  area  where  this  element  is  most  likely  to  be 
useful.  The  number  of  published  problems  available  for  comparison  is  limited,  especially 
in  the  case  of  dynamic  plastic  deformation.  The  solutions  provided  by  this  element  in  the 
plastic  region  need  to  be  verified  by  experimental  data,  and  the  constants  of  Gurson’s  void 
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model  need  to  be  determined  for  a  variety  of  materials  and  porosity  conditions.  In  all,  there 
are  six  material  property  constants  that  must  be  determined  before  an  accurate  solution  can 
be  obtained:  qx,  q2,  qi,f$,  and  %.  These  are  not  readily  available,  but  play  a  vital  role  in 
the  constitutive  equations  of  this  formulation. 

Implementing  this  shell  element  into  a  general  purpose  finite  element  program  will 
allow  easier  comparisons  with  other  element  formulations,  including  solid  brick  elements, 
and  will  facilitate  the  analysis  of  more  complex  problems.  This  work  has  already  begun  for 
a  version  of  DYNA3D,  but  is  not  yet  complete. 

The  method  of  hourglass  mode  control  used  implemented  a  full  version  of  the  same 
formulation  used  for  the  internal  force  calculation.  Efficiency  may  be  increased  by  using  a 
simpler  formulation  to  calculate  the  internal  hourglass  forces.  Along  the  same  lines,  no 
attempt  has  been  made  to  improve  the  efficiency  of  the  code  presented  here  so  that 
readability  may  be  preserved. 
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APPENDIX  A.  PREPROCESSOR  IN  MATLAB  FOR  FEA, 


%  feapre.m  -  Preprocessor  for  fea  program 
%  Program  Limitations : 

%  No  input  of  Initial  Conditions 

%  Boundary  Conditions  must  be  applied  manually  to  result  file 

%  Load  must  be  applied  manually  to  result  file 

%  Structure  will  be  homogenous 

%  Only  1/4  cylinder  supported  at  this  time 

%  Get  Input  data 

type  =  input {'Type  of  Structure  (l=cylinder,  2=plate,  3=shoe  box, 
4-capl,  5=cap2) :  ' ) ; 

%  Structure  Specific  Input 
if  type==l 

nelx  =  input (' Number  of  Elements  along  x-axis:'); 

nelp  =  input ( 'Number  of  Elements  perpendicular  to  x-axis:'); 

xl  =  input ('Lower  X-axis  value:'); 

x2  =  input ( 'Upper  X-axis  value:'); 

R  =  input ( ' Radius : ' ) ; 

thetal  =  input ( 'Starting  Angle  (degrees):'); 
theta2  =  input (' Ending  Angle  (degress):'); 

elseif  type==2 

nelx  =  input ( 'Number  of  Elements  along  x-axis:'); 

nelp  =  input ( 'Number  of  Elements  along  y-axis:'); 

xl  =  input ('Lower  X-axis  value:'); 
x2  =  input ('Upper  X-axis  value:'); 
yl  =  input ('Lower  Y-axis  value:'); 
y2  =  input ('Upper  Y-axis  value:'); 

elseif  type==3 

nelx  =  input ( 'Number  of  Elements  along  x-axis:'); 

nelp  =  input ( 'Number  of  Elements  Perpendicular  to  x-axis:'); 

xl  -  input ( ' Lower  X-axis  value : ' ) ; 
x2  =  input ('Upper  X-axis  value:'); 
yl  =  input ('Lower  Y-axis  value:'); 
y2  =  input ('Upper  Y-axis  value:'); 
zl  =  input ('Lower  Z-axis  value:'); 
z2  =  input ('Upper  Z-axis  value:'); 

elseif  type==4 

R  =  input ( ' Radius : ' ) ; 

alpha  =  input ('Phi  (angle  from  Z-axis  of  inscribed  circle):'); 
neap  -  input {'Size  (1=25  elem  2=42  elem  3=63  elem) : ' ) ; 
if  ncap==l 
narc=6; 
ntheta=ll; 
elseif  ncap==2 
narc=8 ; 
ntheta=13; 
elseif  ncap==3 
narc=10; 
ntheta=15; 
else 

errordlg ( ' Bad  Type ' ) ; 

end 
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elseif  type==5 

R  =  input ( ' Radius : ’ ) ; 

alpha 1  =  input ('Alphal  (elevation  angle  to  bottom  ring  (in 
degrees) : ' ) ; 

alpha2  =  input ( 'Alpha2  (elevation  angle  to  top  ring  (in  degrees)):') 
thetal  =  input ( 'Thetal  (start  angle  from  x-axis  (in  degrees)):'); 
theta2  =  input (' Theta2  (end  angle  from  x-axis  (in  degrees)):'); 
nelp  =  input (' Number  of  Elements  along  elevation  angle:'); 
nelx  =  input (' Number  of  Elements  along  theta  angle:'); 
else 

errordlg ( 'Bad  Type'); 

end 

%  Input  Material  and  Time  Data 

E  =  input ( 'Elastic  Modulus:'); 

rho  =  input ( 'Density  (mass  /  unit  volume);!); 

pois  =  input ('Poissons  Ratio:'); 

t  =  input ( ' Thickness : ’ ) ; 

Syp  =  input ('Yield  Stress:’); 

ipt  =  input ( 'New(2)  or  Old(O)  shell  element  (enter  2  or  0) : ' ) ; 
stime  =  input ('Start  Time:’); 
etime  =  input ( ' Stop  Time : ' )  ; 

delt  =  input ( 'Calculation  Step  Time  (delta):'); 
ptime  =  input ('Print  Step  Time:'); 

%  If  old  shell  element,  adjust  density  to  mass  /  unit  area 
if  ipt==0 

rho  =  rho  *  t; 
ipoint  =0; 
else 

ipoint  =  input ( 'Number  of  Integration  Points  in  Z  direction:'); 

end 

%  ***  Cylinder  *** 
if  type==l 

%  Calculate  Node  x, y, z  coordinates 
thetal  =  thetal  *  pi  /  180; 
theta2  =  theta2  *  pi  /  180; 
nndsx  =  nelx+1; 
nndsr  =  nelp+1; 
xstep  =  (x2  -  xl)  /  nelx; 
rstep  =  (theta2  -  thetal)  /  nelp; 
i  =  1; 

ne  =  nelx  *  nelp; 
nnodes  =  nndsx  *  nndsr; 

%  Initialize  Vectors 
x  =  zeros (nnodes, 1) ; 
y  =  x; 
z  =  x; 

nodes  =  zeros (ne, 4); 

for  xt=xl: xstep :x2 

for  r=thetal: rstep :theta2 
x ( i )  =  xt ; 
y (i)  =  -R  *  cos (r) ; 
z (i)  =  R  *  sin (r) ; 
i  =  i  +  1; 

end 

end 
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xd  =  x  /  R; 
yd  =  y  /  R; 
zd  =  2  /  R; 

end 

%  ***  Plate  *** 
if  type==2 

nndsx  =  nelx+1; 
nndsy  =  nelp+1; 
xstep  -  (x2  -  xl)/nelx; 
ystep  =  (y2  -  yl)/nelp; 
i  «  1; 

ne  =  nelx  *  nelp; 
nnodes  =  nndsx  *  nndsy; 
x  =  zeros (nnodes, 1) ; 
y  =  x; 
z  =  x; 

nodes  =  zeros (ne, 4); 
xd  =  zeros (nnodes, 1) ; 
yd  =  xd; 

zd  -  -I  *  ones (nnodes, 1) ; 

for  xt  =  xl: xstep :x2 

for  yt=yl: ystep :y2 
x(i)  =  xt; 
y(i)  »  yt; 
i  =  i  +  1; 

end 

end 

end 

%  ***  Shoe  Box  *** 
if  type=-3 

nndsx  =  nelx  +  1; 
nndsp  =  nelp  +1; 
xstep  «  (x2  -  xl)/nelx; 
ystep  =  (y2  -  yl)/(nelp  /  2); 
zstep  =  (z2  -  zl)/(nelp  /  4); 
i  =  1; 

ne  =  nelx  *  nelp  ; 
nnodes  =  nndsx  *  nndsp; 
x  =  zeros (nnodes, 1) ;  ' 
y  =  x; 
z  =  x; 
xd  =  x; 
yd  =  x; 
zd  =  x; 

nodes  =  zeros (ne, 4); 

for  xt  =  xl: xstep :x2 
%  Front 

for  zt  =  zl: zstep: z2 
x (i)  =  xt; 
y  (i)  =  yi; 

yd(i)  -  1; 

z  ( i )  -  zt ; 
i  =  i  +  1; 

end 
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%  Top 

for  yt  =  yl+ystep:ystep:y2 
x  (i)  =  xt; 
y  (i)  -  yt; 
z(i)  =  z2; 
zd(i)  =  -I; 

i  =  i  +  1;  • 

end 
%  Back 

for  zt  -  z2-zstep:-zstep: zl 
x  (i)  =  xt; 
y  (i)  =  y2; 
yd(i)'  -  -1; 

z(i)  =  zt; 
i  =  i  +  1; 

end 

end 

end 

%  Spherical  Cap  #1 
if  type==4 

%  Generates  FEA  mesh  of  spherical  cap  using  symmetry 
%  Also  handles  nodal  connectivity 

%  Input  Values 

theta  =  linspace (0, pi/2, ntheta) ; 
alpha  =  pi  *  alpha  /  180; 

phi  =  linspace (0, alpha ,narc) ; 

%  Spherical  Coordinates  for  bottom  nodes 
for  i=l:ntheta 

P (i, : )  =  [R  theta (i)  phi(narc)]; 
bct(i)  =  7; 
bcr ( i )  =  7  ; 
end 

%  Spherical  Coordinates  for  next  row  of  nodes 

nnodes  =  ntheta; 

for  i  -  1 : (ntheta-1) /2 

P (i+ntheta, : )  =  [Rtheta(2  *  i)  (phi (narc) - (phi (narc) 
phi (narc-1) ) /2) ] ; 

end 

%  Spherical  Coordinates  for  the  rest  of  the  nodes 

nnodes  =  nnodes  +  (ntheta-1 ) /2 ; 

nrow  =  ntheta  -  (ntheta-1) /2; 

for  j=narc-l:-l: 1 

Th  =  0; 

dTh  =  pi  /  (2  *  (nrow  -  1) ) ; 
bet (nnodes+l)=2; 
bcr (nnodes+l)=6; 
for  i=l : nrow 

P (i+nnodes, : )  =  [R  Th  phi(j)]; 

Th  =  Th  +  dTh; 

end 

nnodes  =  nnodes  +  nrow; 
bet (nnodes) =1; 
bcr (nnodes) =5; 

nrow  =  nrow  -  1; 

end 
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%  Change  BC  of  top  node 
bet (nnodes) =4; 
ber (nnodes) =7; 

%  Convert  to  Cartesian  Coordinates 
for  i=l: nnodes 

x(i)  =  P(i,l)  *  cos (P (i, 2) )  *  sin (P (i, 3) ) ; 
y (i)  =  P(i,l)  *  sin (P (i, 2) )  *  sin(P(i,3)); 
z(i)  =  P (i/ 1)  *  cos (P (i, 3)  ) ; 

end 

%  Element  Connections 
%  Bottom  row  first 
mid  =  ntheta; 

top  =  mid  +  ntheta  -  (ntheta  -  l)/2  -1; 
ne  =  1; 

for  i=l: (ntheta-1) /2, 
cl  =  (i— 1 ) *2  +1; 

nodes (ne, : )  =  [cl  (cl+1)  (mid+i)  (top+i)]; 
nodes (ne+1, : )  =  [(cl+1)  (cl+2)  (top+l+i)  (mid+i)] 
ne  =  ne  +  2; 

end 

%  Rest  of  elements 

nrow  =  ntheta  -  (ntheta-1) /2  -  1; 

bot  =  mid; 

mid  -  top; 

top  =  top  +  nrow  +  1; 
for  j=l:nrow- 
for  1=1: nrow 
cl  =  mid  +  i; 

nodes (nef:)  =  [cl  (bot+i)  (cl+1)  (top+i)]; 

ne=ne+l; 

end 

bot=mid+l; 
mid=top; 
t  op= t  op+nr ow ; 
nrow=nrow-l; 

end 

ne  =  ne  -  1; 

end 


%  Spherical  Cap  #2  (Has  Hole  it  top) 
if  type~-5 

nndsx  '=  nelx+1; 
nndsy  =  nelp+1; 
theta2  =  theta2  *  pi  /  180; 
thetal  =  thetal  *  pi  /  180; 
alpha2  =  alpha2  *  pi  /  180; 
alphal  =  alphal  *  pi  /  180; 

tstep  =  (theta2  -  thetal) /nelx; 
astep  =  (alpha2  -  alphal) /nelp; 
i  =  1; 

ne  =  nelx  *  nelp; 
nnodes  =  nndsx  *  nndsy; 
x  =  zeros (nnodes, 1) ; 
y  =  x; 
z  =  x; 
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nodes  -  zeros (ne, 4); 


for  thetat  =  thetal : tstep: theta2 

for  alphat=alphal:astep:alpha2 

x (i)  =  R  *  cos (alphat )  *  cos (thetat); 
y (i)  =  R  *  cos (alphat)  *  sin (thetat); 
z  (i)  =  R  *  sin (alphat); 
i  =  i  +1; 

end 

end 

end 

%  Calculate  Element  Connectivity  (will  work  for  all  shapes  except  cap) 
if  type~=4 

for  i=0: (nelx  -  1) 

for  j=l:nelp 

ie  -  (i  *  nelp)  +  j; 
nodes(ie,l)  =  (i  *  (nelp  +1))  +  j; 
nodes (ie, 2)  =  nodes (ie,l)  +  nelp  +  1; 
nodes (ie, 3)  =  nodes (ie, 2)  +  1; 
nodes (ie, 4)  =  nodes (ie,l)  +  1; 

end 

end 

end 

%  Void  content 

PhiO  =  input ( ’Initial  Void  Content: 1 ) ; 
ql  =  input ('Value  for  ql : ' ) ; 

q2  =  input ('Value  for  q2  (Enter  0  for  Default):'); 
q3  =  input ('Value  for  q3  (Enter  0  for  Default):'); 
fn  =  input ( 'Nucleating  Particle  Content:'); 

en  =  input ('Mean  Nucleation  Strain  (Enter  0  for  Default) :'); 

sn  =  input ( ’Nucleation  Standard  Deviation  (Enter  0  for  Default):'); 

%  Create  Stress-Strain  curve 
correct  =  'N'; 

while  ((correct  —  'n')  I  (correct  ==  'N')) 

iseg  =  input ( 'Number  of  Plastic  Segments  (0-9) :') ; 
if  (iseg  -=  0) 
for  i=l:iseg 

slope (i) =input (' Slope  of  Section:'); 

strain ( i ) -input ( ' Right  Sr rain  limit  of  Section:'); 

end 

end 

strn(l)  =  0; 
strss{l)  =  0; 
strn (2)  =  Syp/E; 
strss(2)  =  Syp; 
if  (iseg~=0) 
for  i=l:iseg 

strn (2+i) =strain (i) ; 

strss (2+i) =strss (2+i-l)  +  (strn (2+i) -strn (2+i-l) ) *slope (i) ; 

end 

end 

figure (1) ; 

plot (strn, strss) ,grid,xlabel ( 'Strain' ) , ylabel ( 'Stress’ ) , title ( 'Mat 
erial  Yield  Property') 

correct  =  input ('Is  the  Stress-Strain  Plot  Correct?  (y  or  n):','s’); 
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end 


%  Fill  eprop 
eprop=zeros  (1,40); 
eprop (1)  =  rho; 
eprop ( 2 )  =  E ; 
eprop (3)  =  pois; 
eprop ( 4 )  =  t ; 
eprop (5)  =  Syp; 
eprop (6)  =  ql; 
eprop  (7)  =  q2 ; 
eprop (8)  =  q3; 
eprop (9)  =  fn; 
eprop (10)  -  en; 
eprop (11)  -  sn; 
eprop (12)  =  PhiO; 
eprop (13)  =  iseg; 
if  (iseg  >  0) 
for  i=l:iseg 

eprop (12+2*i)=slope(i) ; 
eprop ( 12+2*i+l ) =strain ( i ) ; 

end 

end 

%  Display  resulting  structure 
figure (2) 

AdjView  =  [-37.5  30]; 

shownode; 

rotate3d  on; 

%  Apply  Boundary  Conditions 
if  type~=4 

bet  =  zeros (1, nnodes )  ; 
ber  =  bet; 

end 

disp ( 1  Boundary  Conditions : 1 ) ; 
dispt'O  =  No  Constraint1); 
disp(Tl  =  Constrained  in  X1); 
disp(*2  =  Constrained  in  Y'); 
disp ( 1 3  =  Constrained  in  2 1 ) ; 
disp('4  =  Constrained  in  X  and  Y’); 
disp ( 1 5  =  Constrained  in  Y  and  Z 1 ) ; 
disp(?6  =  Constrained  in  X  and  Z’); 
disp ('7  =  Constrained  in  X,Y,  and  Zf); 
disp ( 1 ....')  ; 
if  type==4, 

disp (f (It  is  inadvisable  to  change  BC,fs  for  Spherical  Cap)') 

end 

nod  =  input (’Node  to  change  Boundary  Conditions  (0  to  end,  -1  to 
list)  :  ’ )  ; 

while  ((nod  0)  &  (nod  <=  nnodes) ) 
if  nod  <  0 

disp ( ’ Node  Trans  Rot 1 ) ; 
for  i=l: nnodes 

disp (sprintf ( ’  %d  %d  %d* ,i,bct (i) ,bcr (i) ) ) ; 

end 

else 

disp  (sprintf  ( ’Current :  Translation  =  %d,  Rotation  = 
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' %d' ,bct (nod) ,bcr (nod) ) ) ; 

bet (nod)  =  input ('New  Translation  Boundary  Condition (0-7 ):') ; 
bcr(nod)  =  input ('New  Rotational  Boundary  Condition (0-7) :') ; 
dispC  '); 

end 

nod  =  input ('Node  to  change  Boundary  Condition  (0  to  end):'); 

end 

%  Create  Load  History 

disp( ' (Pressure  Load  only  works  for  Flat  Plate)'}; 

nload  =  input (' Number  of  nodes  to  apply  load  to  (0  =  pressure  load):') 

ntime  =  input (' Number  of  time  intervals  for  Load  Description: ' ) ; 

disp(' Force  Direction:'); 

disp('l  =  Force  along  X-axis'); 

disp('2  =  Force  along  Y-axis'); 

disp('3  =  Force  along  Z-axis'); 

disp('4  =  Moment  about  X-axis'); 

disp('5  =  Moment  about  Y-axis'); 

dispC  6  =  Moment  about  Z-axis'); 

if  nload  ==  0 

nload  =  nnodes  *  6; 

disp(' Input  Pressure  Time  History'); 
for  i=l: ntime 

ldtime (i)  =  input ( ' Time : ' ) ; 

Pressure  (i)  =  input  ( 'Magnitude  (use  neg.  values  for  bottom',... 
'pressure) : ' ) ; 

end 

F  =  zeros (nnodes, 1) ; 
for  elem=l:ne 

%  Get  Coordinates  of  Nodes  for  this  element 

xl  =  x (nodes (elem, 1) ) ; 

x2  =  x (nodes (elem, 2) ) ; 

x3  =  x (nodes (elem, 3) ) ; 

x4  =  x (nodes (elem, 4) ) ; 

yl  =  y (nodes (elem, 1) ) ; 

y2  =  y (nodes (elem, 2) ) ; 

y3  =  y (nodes (elem, 3) ) ; 

y4  =  y (nodes (elem, 4) ) ; 

zl  =  z (nodes (elem, 1) ) ; 

z2  =  z (nodes (elem, 2) ) ; 

z3  =  z (nodes (elem, 3) ) ; 

z4  =  z (nodes (elem, 4) ) ; 

%  Calculate  length  of  sides,  then  area 
LI  =  sqrt ( (x2-xl) A2+ (y2-yl) A2+ (z2-zl) A2) ; 

L2  =  sqrt ( (x3-x2) A2+ (y3-y2) A2+ (z3-z2) A2) ; 

L3  =  sqrt ( (x4-x3) A2+ (y4-y3) A2+ (z4-z3) A2) ; 

L4  =  sqrt ( (xl-x4) A2+ (yl-y4) A2+ (zl-z4) A2) ; 

dotl  =  (Xl-x2)*(x3-x2)+(yl-y2)*(y3-y2)+(zl-z2)*(z3-z2); 

dot2  =  (xl-x4) * (x3-x4)+(yl-y4) * (y3-y4) + (zl-z4) * (z3-z4) ; 

Area  =  0.5  *  (sqrt(LlA2  *  L2A2  -  dotlA2)  +  sqrt  (L3A2  *  ... 

L4A2  -  dot2A2) ) ; 

%  Calculate  Force  per  unit  pressure  (assumes  square  elements) 
F (nodes (elem, 1) )  =  F(nodes (elem, 1) )  +  Area  /  4; 

F (nodes (elem, 2) )  =  F(nodes (elem, 2) )  +  Area  /  4; 

F (nodes (elem,  3) )  =  F (nodes (elem, 3) )  +  Area  /  4; 

F (nodes (elem, 4 ) )  =  F (nodes (elem, 4 ) )  +  Area  /  4; 

end 
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%  Create  Node, Force, Time  Vectors 

nc  =  1; 

for  n=l:nnodes, 

loadnode(nc) ’=  n; 

loadnode (nc+1)  =  n; 

loadnode (nc+2)  =  n; 

loadnode (nc+3)  =  n; 

loadnode (nc+4 )  =  n; 

loadnode (nc+5)  =  n; 

loaddir(nc)  =  1; 

loaddir (nc+1)  =  2; 

loaddir (nc+2)  =  3; 

loaddir (nc+3)  =4; 

loaddir (nc+4)  =  5; 

loaddir (nc+5)  -  6; 

loadtime (nc, : )  «  ldtime; 

loadtime (nc+1, : )  =  ldtime; 

loadtime (nc+2, : )  =  ldtime; 

loadtime (nc+3, : )  =  ldtime; 

loadtime (nc+4, : )  =  ldtime; 

loadtime (nc+5, : )  -  ldtime; 

loadmag (nc, : )  =  xd(n)  *  F(n)  *  Pressure; 

loadmag (nc+1, : )  =  yd(n)  *  F(n)  *  Pressure; 

loadmag (nc+2, : )  =  zd(n)  *  F(n)  *  Pressure; 

loadmag (nc+3, : )  =  loadmag (nc, : )  *  t  /  2; 

loadmag (nc+4, : )  =  loadmag (nc+1, : )  *  t  /  2; 

loadmag (nc+5, : )  =  loadmag (nc+2, : )  *  t  /  2; 

nc  =  nc  +  6; 

end 

else 

disp ( 1  Input  Load/Time/Magnitude  History* ) ; 

for  n=l:nload 

loadnode (n)  =  input (’Node  Number:’); 

loaddir (n)  =  input (’Load  Direction:’); 

for  nt  =  l:ntime 

loadtime (n,nt)  =  input ( ’Time: ’) ; 
loadmag (n, nt )  =  input ( ’ Magnitude : ’ ) ; 

end 

end 

end 

%  Create  draft  input  file 
F  =  f open ( ’input .in’ , ’w’ ) ; 

fprintf (F, ’ %d  1  %d  %d  0  %d  0  \r\n’ , nnodes, nload, ntime, ipt ) 
fprintf (F, *  %e  %e  %e  %e  \r\n’ , stime, etime, delt,ptime) ; 
fprintf (F, '0  00000  \r\n’); 
fprintf (F, f0  %d  0  0  0  \r\n’,ne); 
fprintf (F, ’0  000  \r\n’); 

'fprintf (F, ’ 0  0\r\n0  0\r\n%d\r\nl  1  \r\n’ , ipoint ) ; 
for  i=l:4 

fprintf  (F,  ’  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e\r\n’,... 

eprop ( (1+ (i-1) *10) : (i*10)  ) )  ; 

end 

for  i=l: nnodes 

fprintf  (F,  ’  %d  %d  %d  %e  %e  %e  \r\n’,... 

i,bct (i) ,bcr (i) ,x(i) , y  (i) , z  (i) ) ; 


for  i=l:ne 

fprintf (F, ' %d  1  %d  %d  %d  %d  0  \r\n' , i, nodes (i, :)) ; 

end 

for  i=l:nload 

fprintf (F, ' %d  %d\r\n' , loadnode (i) , loaddir (i) ) ; 
for  nt  =  l:ntime 

fprintf (F, '%e  %e\r\n' , loadtime (i, nt) , loadmag (i,nt) ) 

end 

end 

fclose (F) ; 
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APPENDIX  B.  POST-PROCESSOR  IN  MATLAB  FOR  FEA. 


%  fea_disp.m 
%  P .  McDermott 
% 

%  Uses  GUI  to  read  and  display  output  from  fea  program 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Scrn  =  get (0, ’ScreenSize’ ) ; 

ht  =  round (0.9  *  Scrn(4)  /  2  -  50); 

wd  =  round (0.9  *  (Scrn (3)  -  120)  /  2); 

hf  =  figure ( ’NumberTitle’ , 1  off ’, ’Name1, f  FEA1 , ’MenuBar',  ’none',  . . . 

’Position* , [20  (Scrn(4)  -  500)  120  460] , ’Color *, ’bf ) ; 

hfs  =  figure ( ’NumberTitle ’ , ’off ’ ,  ’Name’ ,  ’Structure 
View’ , ’Visib’ , ’off’ , . . . 

*  Position’ , [160  (Scrn (4 ) -70-ht)  wd  ht] ) ; 

hfd  =  figure ( ’NumberTitle ’ , ’off' , ’Name ’ , ’ Displacement 1 , ’Visib1 , ’off 
’ Position’ ,[ (wd+180)  (Scrn (4)  -  70-ht)  wd  ht] ) ; 


hlast  =  0; 

AdjView  =  [-37.5  30]; 

flagd  =  1; 

flagv  =  0; 

flaga  =  0; 

flags  -  0; 

File  =  ’output . fea’ ; 

File2  =  ’ output . fes ’ ; 

P  =  pwd; 

Path  =  strcat (P, ’ \ ’ ) ; 

hc_change  =  uicontrol (hf, ’Style’ , ’push’ , ’String' , ’Change  File’,... 
’ Position’ , [20  380  80  20],... 

’CallBack’ , [. . . 

’ [Filel  Pathl]  =  uiget file (’’*. fea Load  FEA  Output 
File  ”);',... 

’if  Filel  ~=  0, ’ . . . 

'File  =  Filel; ' . . . 

'Path  =  Pathl; ’ . . . 

’ File2  =  File; ’ . . . 

’File2  (length  (File2) )  =  ”s”  ' 

'set (hc_file, ' ’String’ ’ ,File) ; ' . . . 

’end’]); 

hc_load  =  uicontrol (hf, ’Style’ , ’check’ , ’String’ , 'Load  File’,... 

*  Position’ ,  [20  420  80  20],... 

’CallBack’, [.. . 

’ FilePath  -  [Path  File];’,... 

’ FilePath2  =  [Path  File2] 

’parse; ’ , . . . 

'set  (hc__clear,  ’  ’Enable’  ’  ’on'  ’);',... 

' set (hc_disp, ' 1  Enable on 1 ');’,.. . 

’set (hcjvelo, ' ’Enable’ ’,’ ’on’ ’);’,.. . 

’set (hc_acc, ' ’Enable’ ’,’’ on’ ’);’,.. . 

’set (hc_show, ’ 'Enable* ’,’’ on’ ’);’,.. . 
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' set (hc_load, ' 'Value' ’,1, ' 'Enable' ', ' 'off' ');']); 

hc_clear  =  uicontrol (hf, 'Style' , 'push' , 'String' , 'Clear 
Data' , 'Enable ' , 'off' ,  . • • 

' Position' , [20  60  80  20] , 'CallBack' , [ . . . 

'clear  n*  x*  y*  z*  i*  s*  t*  v*  a*  fp  factor  d  etime',... 

'dt  m  time; ' , . • . 

'  set  (hfd,  '  ’Visib",  '  'off' 

'if  (flagv>0) , set (hfv, ' 'Visib' ' , ' 'off' ' ) ;end; ’ , . . . 

'if  (flaga>0) , set (hfa, ' 'Visib' ' , ' 'off' ' ) ;end; ’ , . . . 

'  set  (hfs,  » 'Visib",  '  'off' 

' f lagd  =  1;  flagv  =  0;  flaga  =  0;  flags  =  0;',... 

' set (hc_load, ' ' Value  "  , 0 , ' ' Enable  "  ,  "  on  . 

' set (hc_disp, ' ' Enable ' ' , ' ' off 
' set (hc_velo, ' ' Enable ' ' , ' ' off 
' set (hc_acc, ' ' Enable ' ' , ' ' of f — 

'  set  (hc_show,  '  '  Enable '  ' ,  '  '  of  f  . 

' set  (hc_replay,  '  'Enable' ' ,  '  'off'  ');',... 

'hlast  =  0;  ' ,  .  . . 

' set (hc_mview, ' ' Enable ' ' , ' ' off 
' set (hc_clear, ' ' Enable  "  , ' ' off  ");']); 

hc_close  =  uicontrol (hf, 'Style', 'push', 'String', 'Close', . .. 

'Callback' , [ . . . 

'close (hfd) ; ' , — 

'if  (flaga>0) , close (hfa) ;end; ' , . . . 

' if  (flagv>0) , close (hfv) ;end; ' , . . . 

'close (hfs) ; ' , — 

' clear; ' , . . . 

'close; '],... 

' Position' , [20  20  80  20]); 

hc_disp  =  uicontrol (hf, 'Style', 'push',  'String', 'Displace', ’Position',. 
[20  260  80  20] ,  . . . 

'Enable' , 'off' ,  . . . 

'Callback' , [ • . . 

'if  (flags) ' , . . . 

' AdjView  =  get (hax, ' 'View' ');',... 
end;  , « • • 

' figure (hfd) 

'displace; ' ,  . .  . 

' set (hc_replay, ' ' Enable ’ ' , ' ' on  ");',.. . 

' flagd  =  0; ' , . . . 

'hlast  =  hfd; ' ] ) ; 

he  mview  =  uicontrol (hf, ' Style ', 'push' ,' String' ,  'Adjustment',... 

~  ' Position' , [20  100  80  20],... 

'Enable' , 'off' , — 

'CallBack' , [ . . . 

' figure (hfs) ;',... 

'mmview3d; ' ] ) ; 

hc_file  =  uicontrol (hf, 'Style', 'text',  'Position',  [20  340  80  20],... 

' String' , File) ; 
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hc_vel°  =  uicontrol {hf,' Style 1 ,  'push' , f  String1 , '  Graphs  ' ,  '  Position1 ,  [20 
220  80  20] , . . . 

Enable1 ,  *  off 1 ,  . .  . 

'Callback1, [. . . 

1 flagv=l; ' , . . . 

' control; 1 ] ) ; 

hc_acc  =  uicontrol (hf , 'Style' , 'push' ,  'String' , 'Stress' , 'Position* ,  [20 
180  80  20] ,  .  . . 

'Enable' , 'off , . . . 

'Callback' , [ . . . 

'if  (flags) ' ,  .  . . 

' AdjView  =  get (hax, ' 'View' ');',... 

'end; ' ,  . . . 

' figure (hfd) 

' stresscolor; ' , . . . 

' flagd  «  0; 1 , . . . 

' set (hc_replay, ' 'Enable' ' , ' 'on' . 

'hlast  =  hfd; '] ) ; 

hc_replay  -  uicontrol (hf, 'Style', 'push1, 'String', 'Replay', . . . 

'Position' , [20  140  80  20],... 

'Enable', 'off', . . . 

'Callback' , [. . . 

' figure (hlast) ;',... 

'movie  (m) ;  *  ]  )  ; 

hc_show  =  uicontrol (hf, 'Style' , 'push' , 'String' , 'Show' ,.. . 

' Position' , [20  300  80  20],... 

'Enable', 'off', . . . 

'Callback' , [ . . . 

'figure (hfs) 

' shownode; ' , . . . 

'  set  (hc_mview,  '  '  Enable  ",  "on* ');',..  . 

' flags  =  1; ' , . .  . 

'hax  =  gca; ' ] ) ; 

%%%%%%%%%%%%%%%%%%%%%%  End  FEA  DISP.M  %%%%%%%%%%%%%%%%%%%%%%%%%%%%«%%%% 
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%  PARSE. M 

%  This  m-file  reads  the  output  file  generated  by  fea.f 

%  and  parses  it  into  matlab  usable  variables  which  will  be  used 

%  by  another  program  to  graphically  display  the  data  in  various 

%  formats 

%  McDermott  1999 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Create  Progress  Box  and  1st  Message 

hprog  =  figure ('NumberTitle', ’off', 'Name', 'Load  File 

Progress ' , 'MenuBar ' , ' none ' ,  . . . 

'Position*, [ (Scrn (3) /2  -  80)  (Scrn(4)/2  -  40)  260  80] , 'Color' , 'y' ) ; 

hmsg  =  ui control  ( 'Style' , 'text' ,' Position' ,... 

[20  45  220  20] , ' BackgroundColor ' ,  ,'y' ,  . .  . 

’ String ' , ' Opening  Files ...'); 
hmsg2  =  uicontrol  (' Style ',' text ',' Position' , ... 

[20  15  220  20] , 'BackgroundColor ', 'y' . 

' String ' , FilePath) ; 
drawnow; 

fp=f open ( FilePath, ' r ' ) ; 
fp2  =  fopen(FilePath2, 'r' )  ; 

top=fscanf (fp, ' %e' ,  33)  ; 

%  Update  Progress  Box 

set (hmsg, 'String' , 'Loading  Displacement' ) ; 
drawnow; 

%  split  into  var  names 

set (hmsg2, 'String' , 'Input  &  eprop' ) ; 
drawnow; 

nnodes=top(l)  ; 
nmat=top(2); 
nload=top (3) ; 
loadtime=top ( 4 ) ; 

•initc=top  (5) ; 
stime=top (8) ; 
etime=top (9) ; 
dt=top (11) ; 
ne=top(19) ; 
ipoint  =  top (31); 

for  i=l:nmat, 

eprop  =  fscanf (fp, ' %e' , 40) ; 

end; 

for  i=l : (nnodes) , 

fscanf (fp, '%e', 3) ;  %  dummy  index 

x (i) =fscanf (fp, '%e' , 1) ;  %  x  coordinates  of  nodes 

y (i j  =f scanf (fp, ' %e ' , 1) ;  %  y  coordinates  of  nodes 

z (i) =fscanf (fp, ' %e' , 1) ;  %  z  coordinates  of  nodes 

end; 
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%  associate  nodes  w/  elements 


for  i=l:ne, 

fscanf (fp, 1 %e' ,2) ; 
nodes (i, : ) -fscanf (fp, 1 %e 1 ,  4) ' ; 
fscanf (fp, '  %e' , 1) ; 

end; 

%initial  conditions 
fscanf (fp, 1 %ef ,  ( (3* (6*nnodes) )  +  (2*nload* (loadtime+1)  )  )  )  ; 
nst eps=round ( ( et ime-st ime ) /dt ) ; 
time  =  []; 
temp  -  nsteps; 
for  t=l insteps, 

set  (hmsg2,  ’String* ,[ 1  Time  Step:  *  num2str(t)  ?  /  *  num2str  (nsteps)  ])  ; 
drawnow; 

test=f scanf (fp, 1 %e 1 , 1) ; 
if  (test  >0.0), 
time(t)  =  test; 

for  i=l innodes, 

fscanf (fp, ’ %e' , 1) ;  %dummy  index 

xu (i, t ) =f scanf (fp, ' %e ' , 1) ; 
xv (i, t) “fscanf (fp, f%e',l) ; 
xa (i, t) =f scanf (fp, *  %e ' , 1) ; 
fscanf (fp,f%e’,l); 
yu(i,t)=f scanf (fp, ' %e* , 1) ; 
yv(i, t)=f scanf (fp, ' %e’ , 1) ; 
ya (i, t ) -f scanf (fp, 1 %e ' , 1) ; 
fscanf (fp, f%e' , 1) ; 
zu (i, t) =f scanf (fp, *%e’,l); 
zv(i,t)=f scanf (fp, *  %e* , 1) ; 
za (i, t) -f scanf (fp, ’%e’ , 1) ; 
fscanf (fp, '%e\l); 
thxu (i, t) =fscanf (fp, ' %e ' , 1) ; 
thxv (i, t) =f scanf (fp, *  %e 1 , 1) ; 
thxa (i, t) =f scanf (fp, 1 %e 1 , 1)  ; 
fscanf  (fp,  '%e',l); 
thyu (i,t)=f scanf (fp, ’%e',l); 
thyv (i, t) =f scanf (fp, ' %ef ,  1)  ; 
thya (i, t) =f scanf (fp, *  %e 1 , 1)  ; 
fscanf (fp, f%e’,l); 
thzu (i, t) =f scanf (fp, ' %e',l) ; 
thzv  (i,  t )  =f scanf  (fp,'%e',l); 
thza (i, t) =f scanf (fp, ’ %ef , 1)  ; 
end; 
else 

%  EOF  reached,  terminate  loop  storing  maximum  number  of  time  steps 
t  =  nsteps; 
end; 
end; 

nsteps  =  size(zu,2); 
fclose (fp) ; 

%  Update  Progress  box 

set (hmsg, 1  String1 , 'Calculating  Scales  * ) ; 

set (hmsg2, f  String1 , 1  ' )  ; 

drawnow; 

%  Autoscale  vectors 

dx  =  abs(max(x)  -  min(x)); 

dy  =  abs(max(y)  -  min(y)); 
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dz  =  abs(max(z)  -  min(z)); 

%  Check  for  flat  structures 
if  dx  —  0 

dx  =  1; 

end 

if  dy  ==  0 

dy  =  1; 

end 

if  dz  ==  0 

dz  =  1; 

end 

scale (1)  =  dx  /  max (max (abs (xu) ) ) 
scale (2)  -  dy  /  max (max (abs (yu) ) ) 
scale (3)  *  dz  /  max (max (abs (zu) ) ) 

dscale  =  min (scale)  *  0.25; 
if  dscale  <  1 

dscale  =  1; 

end 

xu  =  xu  *  dscale; 

yu  =  yu  *  dscale; 

zu  =  zu  *  dscale; 

scale (1)  =  dx  /  max (max (abs (xv) ) ) 
scale (2)  =  dy  /  max (max (abs (yv) ) ) 
scale (3)  =  dz  /  max (max (abs (zv) ) ) 

vscale  =  min (scale)  *  0.75; 
if  vscale  <  1 

vscale  =  1; 

end 

xv  =  xv  *  vscale; 

yv  =  yv  *  vscale; 

zv  =  zv  *  vscale; 

scaled)  -  dx  /  max  (max  (abs  (xa)  )  ) 
scale (2)  -  dy  /  max (max (abs (ya) ) ) 
scale (3)  -  dz  /  max (max (abs (za) ) ) 

ascale  =  min (scale)  *  0.5; 
if  ascale  <  1 

ascale  =  1; 

end 

xa  =  xa  *  ascale; 

'  ya  =  ya  *  ascale; 

za  =  za  *  ascale; 


%  Allocate  arrays  for  stress-strain  data 

stress  =  zeros (ipoint*ne*6,nsteps) ; 

strain  -  stress; 

plastrn  =  stress; 

void  -  zeros (ipoint*ne, nsteps) ; 

yieldstrs  -  void; 

effplstrn  =  void; 

vmstr ess  =  void; 

effstrain  =  void; 

estrain  =  void; 

%pstrain  =  void; 

%  Update  Progress  Box 

set (hmsg, 'String' , ’  Loading  Stress  &  Strain1); 
drawnow; 

%  Read  Stress-Strain  Data 
for  t=l insteps, 

test  =  fscanf (fp2, ’%s ’,  1) ; 

set (hmsg2, 'String1 ,[ 'Time  Step:  f  num2str(t)  T  /  1  num2str (nsteps) ]) ; 
drawnow; 

if  (test  >  0.0), 

tl=fscanf (fp2, ’%e’,l) ; 
for  elem=l:ne, 

fscanf  (fp2,  *  %s  * , 1) ; 

fscanf  (fp2,  '%dM); 
for  point=l : ipoint 

index  -  (elem-1) *6*ipoint  +  (point-1) *6; 
fscanf (fp2, ’%d’,l) ; 
for  i  =  1: 6, 

ind  =  i  +  index; 

stress (i  +  index, t) =fscanf (fp2, ’ %e ’, 1) ; 
strain(i  +  index, t)=fscanf (fp2, ’%e’ , 1) ; 
plastrn (i  +  index, t ) -f scanf (fp2, ’%e’ , 1) ; 

end 

indexl  =  (elem-1) *ipoint  +  point; 
yieldstrs (indexl, t)  =  fscanf (fp2, ’ %e* , 1) ; 
void (indexl, t)  =  fscanf (fp2, * %e *, 1) ; 
effplstrn ( indexl, t)  =  fscanf (fp2, * %e 1 , 1) ; 

end 

end 

else 

t  =  nsteps; 

end 

end 

%  Calculate  Von-Misses  Stress  and  an  equivalent  total  strain  value 
%  Update  Progress  Box 

set (hmsg, 1  String1 , ’Calculating  Effective1 ) ; 
set (hmsg2, ’String’ , ’Stress  and  Strain...’); 
drawnow; 
f close (fp2) ; 

%  Adjust  for  total  strain 
%totalstrn  =  strain  +  plastrn; 


79 


%  Get  Effective  stress  and  strain  (doesn't  work  with  multiple  materials) 
for  elem=l:ne, 

for  point=l : ipoint , 

index  =  (elem-1) *6*ipoint  +  (point-1) *6  +  1; 
indexl  =  (elem-1) *ipoint  +  point; 

vmstress (indexl, : )=sqrt (0.5* ( (stress (index, : )  -  ... 

stress (index+1, :)) .A2  +  ... 

(  stress  (index+1,  : )  -  stress  (index+2,  :)).  A2  +... 

(  stress (index+2,  : )  -  stress  (index,  :)).  A2  +  6  *... 

(stress  (index+3,  : )  .  A2  +  stress (index+4,  :).  A2  ... 

+  stress (index+5, : ) . A2) ) ) ; 


estrain  (indexl,  :)=sqrt  (2*  ( (strain (index, :) -strain  (index+1,  : ) )  .  A2... 

+ (strain (index+1, :) -strain (index+2) ) .A2  +  (strain  (index+2)... 

-  strain  (index,  :))  .A2)  +  3*  (strain  (index+3)  .  A2  +  ... 

strain ( index+4 ) .A2  +  strain (index+5) . A2) )/ (2* (1+eprop (3) )) ; 

end 

end 

close (hprog) ; 

%%%%%%%%%%%%%%%%%%%%  End  of  PARSE. M  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%  Control. m 

%  Display  Von-Misses  Stress  and  Yield  Stress  curves 
%  McDermott  1999 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


hfv  =  figure ( ’NumberTitle ’ , ’off1 , ’Name’ , 1  Graph  Control ', ’Visib’ ,' on’ . 

’MenuBar’ , 'none1 , 'Color 1 ,  ’b* , ’  Position* , [160  50  wd  ht]  )  ; 

hfa  =  figure ( ’NumberTitle f ,  ’  of f’ , ’Name ’ , 1  Plot  Window’ , ' Visib 1 , 1  of f ’ ,  .  .  . 
1  Position* ,[ (wd+180)  50  wd  ht]  ) ; 

elementstr  -  ['ll1]; 
for  i=2:ne 

elementstr  =  strcat (elementstr, num2str(i) ) ; 
elementstr  =  strcat (elementstr, 1 | 1 ) ; 

end 

pointstr  =  [’ll1]; 
for  i=2 : ipoint 

pointstr  =  strcat (pointstr, num2str (i) ) ; 
pointstr  =  strcat (pointstr,  f | 1 )  ; 

end 

nodestr  =  [ 1 1 1 * ] ; 
for  i=2:nnodes 

nodestr  =  strcat (nodestr, num2str (i)  ) ; 
nodestr  =  strcat (nodestr, 1 | 1 ) ; 

end 

typestr  =  [‘Stress  vs  Strain | Stress  vs  Time [Strain  vs  time  I  Porosity  vs 
Time | 1 , . .  . 

‘Porosity  vs  Strain | Sigma  zz  vs  Thickness | Sigma  xx  vs  Time  I',... 
’Sigma  yy  vs  Time  [Sigma  zz  vs  Time  |  Sigma  Multiplot  | 

’X  Disp.  vs  Time|Y  Disp.  vs  Time  |  Z  Disp.  vs  Timel’,... 

’Theta  x  vs  Timel Theta  y  vs  Timel’,... 

’Theta  z  vs  Time [MeanDisp *  vs  Time | Plastic  Strain | Elastic  Strain’] 
%  Start  at  Bottom 

hvjplot  =  uicontrol  (hfv, ’Style’ , ’push1 , ’Position’,... 

[  (wd/2  -  40)  20  80  20]  ,  1  String’ ,’ Plot 1 ,’ Callback’ , ... 
’flaga=l;drawstress; ’ ) ; 

hv_typetxt  =  uicontrol (hfv, ’Style’ , ’text’ , ’Position’ , [  20  60  120  20],.. 
’String’ , ’Graph  Type:’); 

hv_type  =  uicontrol  (hfv, ’Style’ , ’popupmenu’, ’Position’ ,... 

[  160  60  160  20  ], ’String ’, typestr, ’Value’ ,1,... 

’Callback’ , ’type  =  get (hv_type, ’ ’Value ’’);’) ; 

type  =  1; 

hv_pointtxt  =  uicontrol  (hfv,  ’ Style ’,*  text ',’ Position’ , ... 

[  20  100  120  20  ], ’String’ , ’Integration  Point:’); 

hv_point  =  uicontrol  (hfv, ’Style’, ’popupmenu’, ’Position’,... 

[  280  100  40  20  ], ’String', pointstr, ’Value’,  1,... 

’Callback’ , ’point  =  get (hv_point, ’ ’Value ’’);') ; 

point  =  1; 

hv_elemtxt  =  uicontrol  (hfv,  ’  Style ’,’ text ’,’ Position’ ,  ... 

[  20  140  120  20  ],’ String’ , ’Element  Number:'); 

81 


hv_elem  =  uicontrol  (hfv,  ' Style ', 'popupmenu' ,' Position' ,... 

[280  140  40  20  ], 'String' ,  elementstr, 'Value' ,  1,... 

'Callback' , 'selem  =  get (hv_elem, ' 'Value '');') ; 

selem  =  1; 

hv  nodetxt  =  uicontrol  (hfv, 'Style' , 'text' , 'Position' ,... 

[20  180  120  20  ],' String ',' Node  Number:'); 

hv  node  =  uicontrol  (hfv,  '  Style ' ,  ’ popupmenu ' ,  '  Position ' , ... 

[280  180  40  20] , 'String' ,nodestr, 'Value' ,  1,... 

'Callback' ,' snode  =  get (hv_node, ' 'Value '');') ; 

snode  =1; 

hv_slidtxt  =  uicontrol  (Jifv, 'Style' , 'text' , 'Position' ,  [20  260  120  20],. 
'String' , 'Time  Slider  Value:  '); 

hv_slidtxt2  =  uicontrol  (hfv, '  Style  ’ ,  '  text ' ,  '  Position ' , ... 

[160  260  80  20] , 'String', 'O' ) ; 

hv  slidtxt3  =  uicontrol  (hfv,  '  Style  ’ ,  '  text ' ,  '  Position  ’ , ... 

[240  260  80  20] ,' String' ,  '  Seconds'); 

hv  slid  =  uicontrol  (hfv,  'Style' ,'  slider' ,  'Position'  ,... 

~  [20  220  (wd  -  40)  20] , 'Min' ,1, 'Max', nsteps, 'Value' ,  1,... 

'Callback' , . . . 

['tstp  =  round ( get (hv_slid, ' 'Value' '));',.. . 

'set (hv_slidtxt2, ' 'String' ' ,num2str (time (tstp) ));']); 
tstp  =1; 

set (hfv, ' Visib' ,  'on' ) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%  End  of  CONTROL. M  %%%%%%%%%%%%%%%%%%%%%%%% 
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%  displace. m 

%  Display  the  element  and  its  displacement  at  each  time  step 
%  McDermott  1999 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  plot  the  element (s) 
if  flagd 

zmax=max ( z ) +max (max ( zu) ) + . 1 ; 
ymax=max (y) +max (max(yu) )+. 1; 
xmax=max (x) +max (max (xu) )  + . 1 ; 
zmin=min ( z ) +min (min ( zu) )  - .  1 ; 
ymin=min (y) +min (min (yu) )-.l; 
xmin=min  (x)  +min  (min  (xu)  )  - .  1; 

AxMax  =  [xmin,  xmax,  ymin,  ymax,  zmin,  zmax] ; 

end; 

m=moviein (nsteps) ; 

%  Display  scaling  factor  in  title  bar 

dispstring  =  ['Displacement:  Scale  =  T  num2str (dscale, ' %1 . 4g' ) ] ; 
set(gcf, 'Name' ,dispst ring) ; 

for  t-1: nsteps 

%  display  global  node  numbers 

cla/Xlabel ( ’x  axis' ) rylabel ( fy  axis ' ) , zlabel ( ' Z  axis'), grid  on,, 
set (gca, 'Box' , 'off ') ,hold  on; 

for  i=l:ne 

xl=x (nodes (i, 1) ) +xu (nodes (i,  1) ,  t) ; 
x2=x (nodes (i, 2) ) +xu (nodes (i, 2) , t) ; 
x3=x (nodes (i, 3) ) +xu (nodes (i, 3) , t) ; 
x4=x (nodes (i, 4 ) ) +xu (nodes (i,  4) ,  t) ; 
yl=y (nodes (i, 1) ) +yu (nodes (i,  1)  ,  t)  ; 
y2=y (nodes (i, 2) ) +yu (nodes (i,  2) ,  t) ; 
y3=y (nodes (i, 3) ) +yu (nodes (i, 3) , t) ; 
y4=y (nodes (i, 4) ) +yu (nodes (i, 4) ,  t) ; 
zl~z (nodes (i, 1) ) +zu (nodes (i, 1) , t) ; 
z2=z (nodes (i, 2) ) +zu (nodes (i, 2) , t) ; 
z3=z (nodes (i, 3) ) +zu (nodes (i, 3) , t) ; 
z4=z (nodes (i, 4) ) +zu (nodes (i,  4 ) ,  t) ; 

xvec=[xl,x2,x3/x4/xl]  ; 
yvec— [yl,  y2,  y3,  y4,  yl]  ; 
zvec=[zl, z2, z3, z4, zl] ; 

plot3 (xvec, yvec, zvec) , view (AdjView) , axis (AxMax) ; 

end; 

hold  off; 

m( : , t)=getframe; 

end; 

movie  (m)  ; 

'%%%%%%%%%%%%%%%%%%%%%%%  End  of  DISPLACE. M  %%%%%%%%%%%%%%%%%%%%%%%%% 
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%  shownode.m 

%  Display  the  element 

%  McDermott  1999 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*%%%%%*%%% 

%  Remove  previous  contents 
cl  a;. 

%  plot  the  element (s) 

offset  =  0.1  *  (max (z)  +  max(y)  +  max(x))/3; 

zmax=max ( z) +of f set ; 
ymax=max (y) +offset; 
xmax=max (x) +of f set; 
zmin=min (z) -offset; 
ymin=min (y) -offset; 
xmin=min (x) -offset ; 

AxMax  =  [xmin,xmax,ymin/ymax,zmin,zmax] ; 

xlabel ( ' x  axis ' ) , ylabel ( ’ y  axis ' ) , zlabel ( ' z  axis 
axis (AxMax) , view (Adj View) , grid,  hold  on; 

for  i=l:ne 

xl=x (nodes (i, 1) ) ; 
x2=x(nodes(i,2)); 
x3=x(nodes(i,3)); 
x4=x(nodes(i,4)); 
yl=y (nodes (i, 1) ) ; 
y2=y (nodes (i, 2) ) ; 
y3=y (nodes (i, 3) ) ; 
y4=y (nodes (i, 4) ) ; 
zl=z (nodes (i, 1) ) ; 
z2=z (nodes (i, 2) ) ; 
z3=z (nodes (i,3)); 
z4=z(nodes(i,4) ) ; 

xvec=[xl,x2,x3,x4,xl] ; 
yvec=[yl,y2,y3,y4,yl]; 
zvec=[zl, z2, z3, z4, zl] ; 

plot3 (xvec, yvec, zvec) ; 

end; 

d=offset  *  0.25; 
for  i=l:nnodes 

text (x(i)+d,  y(i)+(2*d) ,z(i)+d,num2str (i) ) ; 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%  End  of  SHOWNODE.M  %%%%%%%%%%%%%%%%%%%%% 
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%  drawstress.m  -  plots  Von-Mises  stress  and  Yield  stress 
%  McDermott  1999 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


figure (hfa) ; 

set  (hfa,  'Visit*1 ,  'on1 )  ; 

flaga  =  1; 


tstp  =  time  vector  index  (not  all  plots  need) 
point  -  selected  integration  point 
selem  =  selected  element 
type  =  Graph  Type  selected: 

1  =  Stress  vs  Strain 

2  =  Stress  vs  Time 

3  =  Strain  vs  Time 

4  =  Porosity  vs  Time 

5  =  Porosity  vs  Plastic  Strain 

6  -  Z  Stress  vs  Thickness  (at  time (tstp)) 

7  =  X  Stress  vs  Time 

8  =  Y  Stress  vs  Time 

9  -  Z  Stress  vs  Time 

10  =  Stress  Multiplot  (all  4  stress  plots  vs  Time) 

11  =  X  Displacement  vs  Time 

12  =  Y  Displacement  vs  Time 

13  =  Z  Displacement  vs  Time 

14  =  Theta  X  vs  Time 

15  =  Theta  Y  vs  Time 

16  =  Theta  Z  vs  Time 

17  =  Mean  Displacement  vs  Time 

18  =  Plastic  Strain  vs  Time 

19  =  Elastic  Strain  vs  Time 


index  =  (selem-1) *6*ipoint  +  (point-1) *6; 
indexl  =  (selem-1) *ipoint  +  point; 


switch  type 

case  1, 

titlestr  =  ['Element  #'  num2str  (selem) , ... 

'  Integration  Point  #'  num2str (point) ] ; 

plot (effstrain (indexl, : ) , vmstress (indexl, : ) , eff strain (indexl, : ) , . . 

yieldstrs (indexl, : ) , ' — ' ) ,grid,xlabel ( 'Effective  Strain' ) 
ylabel  ( '  Effective  Stress 1 ) ,  title  (titlestr) , ... 
legend ('VME  Stress ',' Yield  Stress ' ) 
maxvm  =  max (vmstress (indexl, :) ) 
maxstrn  =  max (effstrain (indexl, : ) ) 


case  2, 

titlestr  =  ['Element  #'  num2str (selem)  '  Integration  Point  # 
num2str (point) ] ; 

plot  (time,  vmstress  (indexl,  : )  ,  time,  yieldstrs  (indexl,  :),... 

'  — ' )  ,grid, xlabel  ( 'Time  (sec)  '),... 
ylabel ( 'Effective  Stress' ) , title (titlestr) , . . . 
legend ('VME  Stress ',' Yield  Stress'); 
maxvm  =  max (vmstress (indexl, :) ) 
maxys  =  max (yieldstrs (indexl,  :)  ) 
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case  3, 

titlestr  =  ['Element  #'  num2str (selem)  '  Integration  Point  #', 
num2str (point) ] ; 

plot (time, ef f strain (indexl, : ) ) , grid, xlabel ( 'Time  (sec) '),... 

ylabel ( 'Effective  Strain' ),  title (titlestr) ; 
maxstrn  =  max (eff strain (indexl, :) ) 

case  4, 

titlestr  =  ['Element  #'  num2str (selem)  '  Integration  Point  #', 
num2str (point) ] ; 

plot  (time,  void  (indexl,  : ) ) ,  grid,  xlabel  ( 'Time  (sec)  '),... 

title (titlestr) ,  ylabel (' Porosity ')  ; 
maxpor  =  max (void (indexl,  :) ) 

case  5, 

titlestr  =  ['Element  #'  num2str (selem)  '  Integration  Point  #', 
num2str (point) ] ; 

plot  (effplstrn  (indexl,  : ) ,  void  (indexl,  : ) )  ,grid,... 
xlabel ( 'Effective  Plastic  Strain'),... 
ylabel ( ' Porosity' ) , title (titlestr) ; 
maxpor  =  max (void (indexl,  : ) ) 

case  6, 

titlestr  =  ['Element  #'  num2str  (selem)  '  Time: 

num2str (time (tstp) )  '  sec']; 

strs  =  zeros (ipoint, 1) ; 
pnts  =  1: ipoint; 
for  p=l : ipoint 

index3  =  (selem-1) *6*ipoint  +  (p-l)*6  +  3; 
strs(p)  =  stress (index3, tstp) ; 

end 

maxstrs  =  max (strs (p) ) 
minstrs  =  min (strs (p) ) 

plot (strs, pnts ), grid, axis  ij , ylabel (' Integration  Point'),- 
xlabel ('Z  Stress'),  title (titlestr) ; 

case  7, 

titlestr  =  ['Element  #’  num2str (selem)  '  Integration  Point  #', 
num2str (point) ] ; 

plot  (time,  stress  (index+1,  : ) ),  grid, xlabel  ( 'Time  (sec)  '),... 

ylabel ( 'X  Stress '), title (titlestr ) ; 
maxstrs  =  max (stress (index+1, :) ) 
minstrs  =  min (stress (index+1, : ) ) 

case  8, 

titlestr  =  ['Element  #'  num2str (selem)  '  Integration  Point  #', 
num2str (point) ] ; 

plot  (time,  stress  (index+2,  : ) ) ,  grid, xlabel  ( 'Time  (sec)  '),... 

ylabel ( ' Y  Stress ' ) , title (titlestr) ; 
maxstrs  =  max (stress (index+2, : ) ) 
minstrs  =  min (stress (index+2, : ) ) 

case  9, 

titlestr  =  ['Element  #'  num2str( selem)  '  Integration  Point  #', 
num2 st r (point) ] ; 

plot  (time,  stress  (index+3,  :) ) ,  grid,  xlabel  (’ Time  (sec)  ') , .... 

ylabel ('Z  Stress '), title (titlestr) ; 
maxstrs  =  max (stress (index+3,  : ) ) 
minstrs  =  min (stress (index+3, : ) ) 
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case  10, 

titlestr  =  [’Element  #’  num2str (selem)  *  Integration  Point  #’, 
num2str (point) ] ; 
subplot (2,2, 1) 

plot  (time,  stress  (index+1,  : )  ) ,  grid,  xlabel  ( ’Time  (sec)  ’),... 

ylabel ( fX  Stress ’ ) , title (titlestr) ; 
subplot (2,2,2) 

plot  (time,  stress  (index+2,  : )  )  ,  grid, xlabel  ( ’Time  (sec)  * ) ,... 

ylabel ( ’ Y  Stress ’ ) , title (titlestr) ; 
subplot (2, 2, 3) 

plot  (time,  stress  (index+3,  : ) ) ,  grid,  xlabel  ( ’Time  (sec)  ’),... 

ylabel ( ’ Z  Stress ’ ) , title (titlestr) ; 
subplot (2,2,4) 

plot  (time,  vmstress  (indexl,  : ) ,  time,  yieldstrs  (indexl,  :),... 

’  — 1 ) ,  grid, xlabel  ( ’Time  (sec)  ’),... 
ylabel  ( ’Effective  Stress’),... 

legend (’VME  Stress ’ ,’ Yield  Stress ’), title (titlestr) ; 
maxstrs_x  -  max (stress (index+1,  : ) ) 
minstrs_x  =  min (stress (index+1, : ) ) 
maxstrs_y  =  max (stress (index+2, : ) ) 
minstrs_y  =  min (stress (index+2, : ) ) 
maxstrs_z  =  max (stress (index+3, : ) ) 
minstrs_z  =  min (stress (index+3, :) ) 

case  11, 

titlestr  =  [’Node  # ’ , num2str (snode) ] ; 

plot  (time,  xu  (snode,  : )  /dscale)  ,  grid,  xlabel  ( ’  Time  (sec)  ’),... 

ylabel {'X  Node  Displacement’ ) , title (titlestr) ; 
maxdx  =  max (xu (snode, :)) /dscale 
mindx  =  min (xu (snode, :)) /dscale 

case  12, 

titlestr  =  [’Node  #’,  num2str  (snode)  ]-; 

plot  (time,  yu  (snode,  : )  /dscale)  ,  grid,  xlabel  ( ’  Time  (sec)  ’),... 

ylabel ( ' Y  Node  Displacement’ ) , title (titlestr) ; 
maxdy  =  max (yu (snode, :)) /dscale 
mindy  =  min (yu (snode, : ) ) /dscale 

case  13, 

titlestr  ='  [’Node  #’, num2str (snode) ] ; 

plot  (time,  2u(snode,  : )  /dscale) , grid, xlabel  ( ’Time  (sec)  ’ )  ,... 

ylabel (’Z  Node  Displacement ’), title (titlestr) ; 
maxd2  =  max (2U (snode, :)) /dscale 
mindz  =  min ( zu (snode, :)) /dscale 

case  14, 

titlestr  =  [’Node  #’, num2str (snode) ] ; 

plot  (time,  thxu (snode,  : )  ) ,  grid,  xlabel  ( ’Time  (sec)  ’),... 

ylabel ( ’Theta  X’ ) , title (titlestr) ; 
maxdz  =  max (thxu ( snode, :) ) 
mindz  =  min (thxu (snode, :) ) 

case  15, 

titlestr  =  [’Node  #’, num2str (snode) ] ; 

plot  (time,  thyu( snode,  : )  )  ,  grid, xlabel  ( ’Time  (sec)  ’),... 

ylabel ( 'Theta  Y’ ), title (titlestr) ; 
maxdz  =  max (thyu (snode, :) ) 
mindz  =  min (thyu (snode, :) ) 
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case  16, 

titlestr  =  ['Node  #', num2str (snode) ] ; 

plot  (time,  thzu( snode,  : ) ) ,  grid,xlabel  ( 'Time  (sec)  '),... 

ylabel ( ' Theta  Z ' ) , title (titlestr) ; 
maxdz  =  max (thzu (snode, :) ) 

'  mindz  =  min (thzu (snode,  :) ) 

case  17, 

titlestr  =  ['Node  #', num2str (snode) ] ; 

du  =  sqrt  (xu  (snode,  :).  A2  +  yu  (snode,  :).  A2  +  ... 

zu (snode, :). A2) /dscale; 
maxdu  =  max(du) 
mindu  =  min(du) 

plot  (time,  du) ,  grid,  xlabel  (' Time  (sec)  '),... 

ylabel ('Mean  Node  Displacement '), title (titlestr); 


case  18, 

titlestr  =  ['Element  #'  num2str (selem)  '  Integration  Point  #',... 
num2str (point) ] ; 

plot  (time,  effplstrn(indexl,  :) ),  grid,  xlabel  (' Time  (sec)  ')  ,... 

ylabel ( 'Effective  Plastic  Strain' ), title (titlestr) ; 
maxplaststrn  =  max (effplstrn (indexl,  : ) ) 

case  19, 

titlestr  =  ['Element  #',  num2str  (selem)  '  Integration  Point  #',.. 
num2str (point) ] ; 

plot (time, estrain (indexl, : ) ), grid, xlabel (' Time  (sec) '),... 

ylabel ( 'Effective  Elastic  Strain' ), title (titlestr) ; 
maxelaststrn  =  max (estrain (indexl, :) ) 

otherwise, 

disp ( ' Invalid  Graph  Type ' ) ; 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%  End  of  DRAWSTRESS..M  %%%%%%%%%%%%%%%%%%%%%%% 
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Figure  23.  Screen  Capture  of  FEA  Preprocessor. 
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APPENDIX  C.  FILE  LISTING  FOR  FEA  SHELL  ELEMENT  SUBROUTINE. 


c*  elsh41.f 

c*  Contains  the  following  subroutines : 

c*  elsh41  -  element  matrices  for  3-D  degenerated  4-node  shell 

q*********************************************************************** 

subroutine  elsh41  (nnds , xcoord, ycoord, zcoord, nesh4 , ndsh4  , matsh4 , 
nmat,matno, eprop,ndsdof , sdisp, svel,  sforce, 
smass, nsdof , ishel4, stssh4,mpoint,  ipoint, itime) 

*Jc  +  **ieieit**  +  ii:*  +  ***ic**ie**'k'k’k'k‘k'k’k'ki''k’k'k'k'k-k-k'k'kic'k'kik'k'k'k'k'k'kit-k-k'k-k'kJc'k'k'k'k-k'k-k'k 

compute  and  assemble  element  matrices  for  shell  elements  with 
4  nodes  (linear  analysis) 

3-D  degenerated  formulation  including  transverse  shear  and  - 
transverse  normal  deformation 
pressure  formulation 

McDermott  1999 

***********************************************************^*** 
common  /fblk/  fail, arrayl (125) 

dimension  xcoord (nnds) , ycoord (nnds) , zcoord(nnds) , 

ndsh4 (nesh4+l, 4 ) ,matsh4 (nesh4+l) , ishel4 (nesh4+l) , 
matno(nmat) ,eprop(nmat, 40) ,ndsdof (nnds) , 
sdisp (nsdof) , svel (nsdof) , sforce (nsdof) ,  smass (nsdof) , 
stssh4  (inpoint  ,nesh4  *8) 

dimension  eforce (24) , Inode (4) , evel (24), slope (10) , dstrnplst (6) , 
emass (24) , gausx (4 ) , lindex (24), strain (10) , corr (2) , 
array2 (15) , estrain(6) ,  estress (6)  ,s(6),A(2,2) ,  Ainv(2, 2) , 
gausy ( 4 ) , weighx ( 4 ) , weighy ( 4 ) , shapel (3) , deri vl (3) , 
shapef (4) ,derivl(2, 4) , derivg(3,4) ,  aj (3,3) ,  ajinv(3,3) , 
bmtx(6,24) , tmp(8) ,xc (4) , yc (4 ) , zc (4) ,estressg(6) , 
derilg (3, 4) , dl C  3 ) , edisp (24) , edispg(24) , rot6t (6,6), 
d2  (3) ,  vl  (3) ,  v2  (3) ,  v3  (3) ,  gausz  (10) , weighz  (10) , 
bmtxt (24,6) , delt (6) , tl (3, 3) , tlinv (3, 3) , B (2) , 
rot6 (6, 6) , rot3 (3, 3) , estrainp (6) , eforceg (24) , 
weighxh(4) ,weighyh(4) ,weighzh(4) ,eforceh(24) , 
gausxh(4) ,gausyh(4) ,gauszh(4) ,iplane(5) 
c 

data  ngausx,ngausy,ndof,n  /l,  1,  6,  4/ 

data  ngausxh,ngausyh,ngauszh,hourl  / 2,  2,  1,  0.05/ 

data  ref,eps  /0.0,1.0e-4  / 

data  delt  /1.0,  1.0,  1.0,  0.0,  0.0,  0.0/ 

data  iplane  /l,  2,  4,  5,  6/ 


c 

c 

c 


neldof=ndof*n 
ngaus  z=ipoint 

ngauss=ngausx*ngausy*ngausz 
ngaus sh=ngausxh*ngausyh*ngauszh 

***  extract  gauss  points  and  weights  for  numerical  integration 


call  gauss 
call  gauss 
call  gauss 
call  gauss 
call  gauss 
call  gauss 


(ngausx, gausx, weighx) 
(ngausy, gausy, weighy) 

( ngaus z, gausz, weighz) 
(ngausxh, gausxh, weighxh) 
(ngausyh, gausyh, weighyh) 
(ngauszh, gauszh, weighzh) 
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o  o 


c 

*** 

c 

c 

£  ■*  ★  -k 

c 


10 

c 

£  *  ** 
c 


20 


c 

0  ★ 
c 


c 

£  *  St  ★ 

c 


loop  for  number  of  elements 
do  370  ie=l,nesh4 
igausp=0 

initialize  element  internal  force  vector 

do  10  i=l,neldof 
eforceh(i)  =  0.0 
eforce(i)  =  0.0 
continue 

extract  material  property  index  and  material  properties 

imt=matsh4  (ie) 
matt yp=matno ( imt ) 
densty=eprop (imt, 1) 
elast  =  eprop(imt, 2) 
poiss  =  eprop{imt,3) 
thick  =  eprop (imt, 4) 

Syp  -  eprop (imt, 5) 
ql  =  eprop (imt, 6) 
q2  =  eprop (imt, 7) 
q3  -  eprop  (imt, 8) 
fn  =  eprop (imt, 9) 
en  =  eprop  (imt, 10} 
sn  -  eprop (imt, 11) 

PhiO  =  eprop (imt, 12) 
iseg  =  eprop (imt, 13)  +  0.1 
ystrain  =  Syp  /  elast 
strain (1)  =  ystrain 

if  (iseg.gt.0)  then 
do  20  i-1, iseg 

slope (i)  =  eprop (imt, 12  +  i  *  2) 
strain (i+1)  =  eprop (imt, 13  +  L  *  2) 
continue 

endif 

bulk  =  elast  /  (3.0  *  (1.0  -  2.0  *  poiss)) 
shear  =  elast  /  (2.0  *  (1.0  +  poiss) ) 

Default  Values 

if  (ql  .ne.  0.0)  then 

if  (q2  .eq.  0.0)  q2  =  1.0 
if  (q3  .eq.  0.0)  q3  =  ql**2 
if  (fn  .ne.  0.0)  then 

if  (sn  .eq.  0.0)  sn  =  0.1 
if  (en  .eq.  0.0)  en  =  0.3 

endif 

endif 

find  dof  index  of  each  element 

do  30  i=l,n 

Inode ( i ) =ndsh4 ( ie , i ) 
do  30  j=l,ndof 

lindex ( (i-1) *ndof+j ) -ndsdof (Inode (i) ) +j-l 
continue 


30 
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c 

q  •k'k'k 

C 


40 

c 

0  ★  ★  * 
c 


50 

c 

0  •*•■*■*■ 

0  ★** 
c 


c 

0  ★  *  ★ 

c 

c 


1 

1 

1 

c 

0  *** 


c 

0  *** 


c 

0  *** 


load  displacement  for  element 

do  40  i=l,neldof 
edispg(i)  =  sdisp (lindex (i) ) 
evel(i)  =  svel (lindex (i) ) 
continue 

extract  element  nodal  coordinate  values 

do  50  i=l,n 
xc ( i ) =xcoord ( Inode ( i ) ) 
yc ( i ) =ycoord ( Inode ( i ) ) 
zc(i)=zcoord (Inode (i) ) 
continue 

compute  direction  vectors 

(vl  &  v2:  tangent ,  v3:  normal  to  surface) 

dl ( 1 ) =xc ( 2 ) -xc ( 1 ) 
dl(2)-yc(2)-yc(l) 
dl (3) =zc (2) -zc (1) 
d2 ( 1 ) =xc ( 4 ) -xc ( 1 ) 
d2 ( 2 ) =yc ( 4 ) -yc ( 1 ) 
d2 (3) =zc (4) -zc (1) 

Calculation  of  Hourglass  Control  Multiplier 

Need  Surface  Area  of  Element 
all=sqrt (dl (1) **2+dl (2) **2+dl (3) **2) 
al4  -  sqrt  (d2 (1) **2+d2 (2) **2+d2 (3) **2) 

al2  =  sqrt((xc(3)-xc(2))**2+(yc(3)-yc(2))**2+(zc(3)-zc(2))**2) 
al3  -  sqrt ( (xc (3) -xc (4 ) ) **2+ (yc (3) -yc (4 ) ) **2+ (zc (3) -zc (4) ) **2) 
adotl  -  (xc (1) -xc (2) ) * (xc (3) -xc(2) )  +  (yc (1) -yc(2) ) * (yc (3) -yc (2) ) 
4-  f  zc  (1 }  -zc  (2\  )  *  f  zc  (3^  -zc  (2)  ) 

adot2  =  (xc  (1)  -xc  (4) )  *  (xc  (3)  -xc(4) )  +  (yc  (1)  -yc(4 )  )  *  (yc  (3)  -yc  (4) ) 

+  (zc (1) -zc (4) ) * (zc (3) -zc  (4) ) 

area  =  0.5* (sqrt (all**2*al2**2-adotl**2) + 

sqrt (al3**2*al4**2-adot2**2) ) 

hour  =  hourl  *  thick**2  /  area 

v3  is  normal  to  plane 

v3  ( 1 )  =dl  (2)  *d2  (3)  -dl  (3)  *d2  (2) 
v3  (2)  =dl  (3)  *d2  (1)  -dl  (1)  *d2  (3) 
v3  (3)  =dl  (1)  *d2  (2)  -dl  (2)  *d2  (1) 
sum=sqrt (v3 (1) **2+v3 (2) **2+v3 (3) **2) 
v3 (1) =v3 (1) /sum 
v3 (2) =v3 (2) /sum 
v3 ( 3 ) =v3 ( 3 ) / sum 

v2  is  unit  vector  in  direction  of  node  1  to  node  4 
v2 ( 1 )  =  d2(l)/al4 
v2 (2)  =  d2(2)/al4 
v2 (3)  =  d2(3)/al4 

vl  is  orthogonal  to  v2,  in  plane  normal  to  v3 
vl  (1)  =v3  (2)  *v2  (3) -v3  (3)  *v2  (2) 
vl  (2)  =v3  (3)  *v2  (1)  -v3  (1)  *v2  (3) 
vl  (3)  =v3  (1)  *v2  (2)  -v3  (2)  *v2  (1) 
sum  =  sqrt (vl (1) **2+vl (2) **2+vl (3) **2) 
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c 

0  *  ★ 


c 

0  ★  ★  ★ 
c 


60 

c 

0  *** 
c 


70 

c 

0  *  ★  ★ 

c 


vl ( 1 ) =-vl ( 1 ) / sum 
vl  (2)=-vl  (2) /sum 
vl  (3)=-vl  (3)  /sum 

Strain  Transfonnation  Matrix 
rot 6 (1, l)=vl (1) **2 
rot6 (1, 2) =vl (2) **2 
rot6 (1, 3) =vl (3) **2 
rot 6 (1, 4)=vl (1) *vl (2) 
rot 6 (1, 5)=vl (2) *vl (3) 
rot 6 (1,6) =vl (1) *vl (3) 
rot  6(2,1) =v2 ( 1 ) *  *2 
rot 6 (2, 2)=v2 (2) **2 
rot6 (2, 3)=v2 (3) **2 
rot  6(2,4) =v2 ( 1 ) *v2 ( 2 ) 
rot  6(2,5) =v2 ( 2 ) *v2 ( 3 ) 
rot 6  (2,  6)  =v2  (1)  *v2  (3) 
rot 6 (3, l)=v3 (1) **2 
rot6(3,2)=v3 (2) **2 
rot 6  (3,  3)  =v3  (3)  **2 
rot6 (3, 4 ) =v3 (1) *v3 (2) 
rot  6(3,5) =v3 ( 2 ) * v3 ( 3 ) 
rot  6(3,6) =v3 ( 1 ) * v3 ( 3 ) 
rot 6  (4, 1)  =2. 0*vl  (1)  *v2  (1) 
rot 6  (4, 2 )  =2. 0*vl  (2)  *v2  (2) 
rot6  (4, 3)  =2 . 0*vl  (3)  *v2  (3) 
rot6(4, 4)=vl (1) *v2 (2)+v2 (1) *vl (2) 
rot  6(4,5) =vl ( 2 ) * v2 ( 3 ) +v2 ( 2 ) *vl ( 3 ) 
rot  6  ( 4 , 6 )  =vl  ( 3 )  *  v2  ( 1 )  +v2  ( 3 )  *  vl  ( 1 ) 
rot 6  (5, 1)  =2 . 0*v2  (1)  *v3  (1) 
rot 6  (5, 2)  =2 . 0*v2  (2)  *v3  (2) 
rot 6  (5, 3)=2 . 0*v2  (3)  *v3  (3) 
rot 6 (5,4) =v2 (1) *v3 (2)+v3 (1) *v2 (2) 
rot  6(5,5) =v2 ( 2 ) * v3 ( 3 ) +v3 ( 2 ) *v2 ( 3 ) 
rot  6(5,6) =v2 ( 3 ) * v3 ( 1 ) +v3 ( 3 ) *v2 ( 1 ) 
rot 6  (6, 1)  =2 . 0*v3  (1)  *vl  (1) 
rot 6  (6,  2)  =2 . 0*v3  (2)  *vl  (2) 
rot 6  (6,  3)  =2 . 0*v3  (3)  *vl  (3) 
rot6 (6, 4  )=v3 (1) *vl (2) +vl (1) *v3 (2) 
rot6 (6, 5) =v3 (2) *vl (3) +vl (2)  *v3  (3) 
rot 6  (6,  6)  =v3  (3)  *vl  (l)+vl  (3-)  *v3  (1) 

Transpose  of  Transformation  Matrix 

do  60  i=l, 6 
do  60  j =1,6 

rot  6t ( i , j ) =r ot  6 ( j , i ) 
continue 

Rotation  tranformation  matrix,  and  its  inverse 

do  70  i=l,3 
tl  (i,  1)  =  vl  ( i ) 
tl  (i, 2)  =  v2(i) 
tl  (i, 3)  =  v3(i) 
continue 

call  inv3x3(tl,tlinv,det) 

Transform  Rotation  displacement  to  local  coordinates 
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80 

c 

£  iHr  ★  ■* 

c 


c 

0  *** 
c 


c 

0  *** 
c 


do  80  i=l,4 
do  8  0  j  =1 , 3 

edisp{ (i-l)*6+j)=edispg( (i-l)*6+j) 
edisp ( (i-1) *6+3+j ) =edispg ( (i-1) *6+4 ) *tlinv ( j , 1) 
+edispg( (i-1) *6+5) *tlinv(j,2) 

!  +edispg ( (i-1) *6+6) *tlinv ( j , 3) 

continue 

Begining  of  Hourglass  control 

if  (hour. gt. 0.0)  then 
do  140  ix=l,ngausxh 
rc=gausxh (ix) 
wx=weighxh (ix) 
do  140  iy=l,ngausyh 
sc=gausyh (iy) 
wy=weighyh (iy) 
do  140  iz=l,ngauszh 
tc-gauszh(iz) 
wz=weighzh ( i z ) 

compute  shape  functions  and  derivatives 

call  shp2d4  (rc, sc, shapef ) 
call  der2d4  (rc, sc, derivl) 
call  shpld2  (tc,shapel) 
call  derld2  (tc,derivl) 

compute  jacobian  matrix  and  its  inverse 

hz=thick*0 .5* (shapel (2) * (1 . 0-ref) -shapel (1) * (1 . 0+ref ) ) 
hzdt=thick*0.5 

a j (1,1) =deri vl (1,1) *xc ( 1 ) +deri vl (1,2) *xc { 2 ) +derivl ( 1 ,  3 ) *xc ( 3 )  + 
derivl ( 1 , 4 ) *xc ( 4 ) +derivl ( 1 , 1 ) *hz* v3 ( 1 ) + 
derivl ( 1 , 2 ) *hz* v3 ( 1 ) +derivl ( 1 , 3 ) *hz  * v3 ( 1 )  + 
derivl (1, 4) *hz*v3(l) 

aj (2, 1) =derivl (2, 1) *xc (1) +derivl (2, 2) *xc (2) +derivl (2, 3) *xc (3)  + 
derivl (2,4) *xc (4 ) +derivl (2,1) *hz*v3 (1) + 
derivl (2,2) *hz*v3 (1) +derivl (2,3) *hz*v3 (1) + 
derivl (2,4) *hz*v3 (1) 
aj (3, l)=hzdt*v3 (1) 

aj (1, 2) =derivl (1,1) *yc (1) +derivl (1,2) *yc (2) +derivl (1,3) *yc (3) + 

derivl (1,4) *yc (4) +derivl (1,1) *hz*v3 (2) + 
derivl ( 1 , 2 ) *hz* v3 ( 2 ) +der ivl ( 1 , 3 ) *hz* v3 ( 2 ) + 
derivl (1,4) *hz*v3 (2) 

aj (2, 2) =derivl (2, 1) *yc (1) +derivl (2, 2) *yc (2) +derivl (2, 3) *yc (3)  + 
derivl (2,4) *yc (4) +derivl (2,1) *hz*v3 (2) + 
derivl (2,2) *hz*v3 (2) +derivl (2,3) *hz*v3 (2) + 
derivl (2, 4) *hz*v3 (2) 
aj (3, 2) =hzdt*v3 (2) 

aj (1, 3) =derivl (1, 1) *zc (1) +derivl (1, 2) *zc (2) +derivl (1, 3) *zc (3) + 
derivl (1, 4 ) *zc (4 ) +derivl (1,1) *hz*v3 (3) + 
derivl (1, 2) *hz*v3 (3) +derivl (1, 3) *hz*v3 (3) + 
derivl (1,4) *hz*v3 (3) 
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c 


aj (2, 3) =derivl (2,1) *zc (1) +derivl (2, 2) *zc (2) +derivl (2, 3) *zc (3) + 

1  derivl (2, 4) *zc (4 ) tderivl (2, 1) *hz*v3 (3) + 

2  derivl (2,2) *hz*v3 (3) tderivl (2,3) *hz*v3 (3) + 

3  derivl (2, 4) *hz*v3 (3) 
aj (3, 3) =hzdt*v3 (3) 

call  inv3x3  (a j , aj inv, det ) 
c 

c  ***  compute  global  derivatives  and  strain-nodal  displacement  matrix 
c 

do  90  i=l,n 

derivg(l, i)=ajinv(l, 1) * derivl (1, i) +ajinv(l, 2) +derivl (2, i) 
deri vg ( 2 , i ) =a j iny (2,1)* derivl ( 1 ,  i ) +a j inv (2,2) *deri vl ( 2 , i ) 
deri vg ( 3 , i ) =a j inv (3,1) *derivl ( 1 ,  i ) +a j inv (3,2) * derivl ( 2 , i ) 
derilg (1, i) =ajinv(l, 3) *hzdt 
derilg (2, i)=ajinv(2, 3) *hzdt 
derilg (3,i)=ajinv(3, 3) *hzdt 
90  continue 

c 

do  100  i=l,6 

do  100  j=l, 24 

bmtx(i, j)=0.0 

100  continue 

c 

do  110  i=l,n 

il=(i-l)*6+l 

i2-il+l 

i3=i2+l 

i4=i3+l 

i5=i4+l 

i6=i5+l 

gkl=derivg (1, i) *hz+shapef (i) *derilg (1, i) 
gk2=derivg (2, i) *hz+shapef (i) *derilg (2, i) 
gk3=derivg (3, i) *hz+shapef (i) *derilg (3, i) 
c 

bmtx ( 1 , il ) =derivg ( 1,  i ) 

bmtx(l,i4)=gkl*(-v2(l) ) 

bmtx ( 1 , i5 ) =gkl * vl ( 1 ) 

bmtx (1, i6) =gkl*v3 (1) 

bmtx(2,i2)=derivg(2,  i) 

bmtx(2,i4)=gk2*(-v2(2)  ) 

bmtx (2, i5) =gk2*vl (2) 

bmtx (2, i6) =gk2*v3 (2) 

bmtx (3, i3)=derivg(3,i) 

bmtx (3, i4 ) =gk3* (-v2 (3) ) 

bmtx (3, i5) =gk3*vl (3) 

bmtx(3,i6)=gk3*v3 (3) 

bmtx ( 4 , i 1 ) =derivg (2,  i ) 

bmtx (4, i2) =derivg (1,  i) 

bmtx (4, i4)=gk2* (-v2 (1) ) +gkl* (-v2 (2) ) 

bmtx (4, ±5) -gk2*vl ( 1 ) +gkl * vl ( 2 ) 

bmtx (4, i6) =gk2*v3 (1) +gkl*v3 (2) 

bmtx ( 5 , i2 ) =deri vg ( 3 , i ) 

bmtx ( 5 , i 3 ) =der i vg ( 2 , i ) 

bmtx ( 5 , i4 ) =gk3* ( -v2 ( 2 ) ) +gk2  * ( -v2 ( 3 ) ) 

bmtx (5, i5) =gk3*vl (2) +gk2*vl (3) 

bmtx (5, i6) =gk3*v3 (2) +gk2*v3 (3) 

bmtx ( 6 , i 1 ) =deri vg ( 3,  i ) 
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110 

c 

120 

c 

c 

0  *** 
c 

c 

c  *** 
c 

c 

c 

130 

c 

0  ★  ★ 
c 


c 

0  *** 
c 

c 

0  **★ 
c 

c 

140 

c 

0  ★  ★ 

c 

c 

0  -k 

c 


bmtx ( 6 , i3 ) =derivg ( 1 ,  i ) 
bmtx(6,i4)=gk3*(-v2(l) ) +gkl* (-v2 (3) ) 
bmtx ( 6 , i5 ) =gk3* vl ( 1 ) +gkl* vl ( 3 ) 
bmtx (6, 16) =gk3*v3 (1) +gkl*v3 (3) 
continue 

do  120  i=l, 6 
do  120  j=l,neldof 

bmtxt ( j , i ) =bmtx ( i , j ) 

weight  =  wx  *  wy  *  wz 
detwt  =  det  *  weight 

calculate  strain 

call  mmult (6, 24, l,bmtx, edisp, estrainp, 0, 1 . 0) 

Transform  strain  to  local  coordinates  system 

call  mmult (6, 6, 1, rot 6, estrainp, estrain, 0,1.0) 

Subtract  any  prior  plastic  strain 
do  130  i=l, 6 

estrain(i)  =  estrain(i)  -  stssh4 (12+i, issts) 

Calculate  stress  using  plane-strain  formulas 

estress(l)  =  elast/ (1 . 0-poiss**2) * (estrain (1)  +  poiss*estrain (2) ) 
estress(2)  -  elast/ (1. 0-poiss**2) * (poiss*estrain (1)  +  estrain (2)) 
estress(3)  =  elast  *  estrain (3) 
estress(4)  =  shear  *  estrain (4) 
estress(5)  =  5.0  /6.0  *  shear  *  estrain(5) 
estress(6)  =  5.0  /6.0  *  shear  *  estrain(6) 

Transform  stress  back  to  global  coordinates 

call  mmult (6,6, l,rot6t, estress, estressg, 0,1.0) 

find  eforceh  for  this  gauss  point  and  sum  to  element  force 

call  mmult (24, 6, 1, bmtxt, estressg, eforceh, 1, detwt) 

continue 

endif 

end  of  hourglass  loop 


numerical  integration  loop  for  membrane 
icount=0 

do  300  ix=l,ngausx 
rc=gausx (ix) 
wx=weighx (ix) 
do  300  iy=l,ngausy 
sc=gausy(iy) 
wy=weighy (iy) 
do  300  iz=l,ngausz 
tc=gausz (iz) 
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wz=weighz (iz) 
c 

c  ***  increase  counter 
c 

igausp=igausp+l 
issts= (ie-1) *8+igausp 
c 

c  ***  Initialize  stssh4  if  first  time 
c 

if  (itime  .  eq.  1)  then 
do  150  i«13,2l' 

150  stssh4 (if issts)  -  0.0 

c 

stssh4 (19, issts)  =  Syp 
stssh4 (20, issts)  -  PhiO 

endif 

c 

c  ***  compute  shape  functions  and  derivatives 
c 

call  shp2d4  (rc, sc, shapef ) 
call  der2d4  (rc, sc,derivl) 
call  shpld2  (tc,shapel) 
call  derld2  (tc,derivl) 
c 

c  ***  compute  jacobian  matrix  and  its  inverse 
c 

hz=thick*0.5* (shapel (2) * (1. 0-ref ) -shapel (1) * (1. 0+ref )  ) 
hzdt=thick*0.5 


a j ( 1 , 1 ) “derivl (1,1) *xc ( 1 ) tderivl (1,2) *xc ( 2 ) +derivl ( 1 , 3 ) *xc ( 3 ) + 

1  derivl ( 1 , 4 ) *xc ( 4 ) +derivl (1,1) *hz*v3 ( 1 ) + 

2  derivl (1, 2) *hz*v3 (1) +derivl (1, 3) *hz*v3 (1) + 

3  derivl (1,4) *hz*v3 (1) 

aj (2, 1) +derivl (2, 1) *xc (1) +derivl (2, 2) *xc (2) +derivl (2,3) *xc(3)+ 

1  derivl (2,4) *xc (4 ) +derivl (2,1) *hz*v3 (1)  + 

2  derivl (2,2) *hz*v3 (1) +derivl (2,3) *hz*v3  (1)  + 

3  derivl (2,4) *hz*v3 (1) 
aj (3, 1) =hzdt*v3 (1) 

a  j  ( 1 , 2 )  =derivl  ( 1 , 1 )  *yc  ( 1 )  tderivl  (1,2)  *yc  (2)  +*derivl  ( 1 , 3 )  *yc  ( 3 )  + 

1  derivl (1,4) *yc (4 ) +derivl (1,1) *hz*v3 (2) + 

2  derivl (1,2) *hz*v3 (2) +derivl (1,3) *hz*v3 (2) + 

3  derivl (1,4) *hz*v3 (2 ) 


aj (2, 2) =derivl (2, 1) *yc (1) +derivl (2, 2) *yc (2) tderivl (2, 3) *yc (3)  + 

1  derivl (2, 4 ) *yc (4 ) +derivl (2, 1) *hz*v3 (2) + 

2  derivl (2, 2) *hz*v3 (2) +derivl (2,3) *hz*v3 (2) + 

3  derivl (2, 4) *hz*v3 (2) 
aj (3, 2) ~hzdt*v3 (2) 

aj (1,3) =derivl (1, 1) *zc (1) tderivl (1, 2) *zc (2) +derivl (1,3) *zc(3)+ 

1  derivl (1, 4) *zc (4 ) tderivl (1,1) *hz*v3 (3) + 

2  derivl (1,2) *hz*v3 (3) +derivl (1, 3) *hz*v3 (3) + 

3  derivl (1,4) *hz*v3 (3) 

aj (2,3) =derivl (2,1) *zc (1) +derivl (2, 2) *zc (2) +derivl (2,3)*zc(3)+. 

derivl (2,4) *zc (4 ) +derivl (2, 1) *hz*v3 (3) + 
derivl (2,2) *hz*v3 (3) +derivl (2,3) *hz*v3 (3) + 
derivl (2,4) *hz*v3 (3) 
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aj (3,3)=hzdt*v3(3) 
c 

call  inv3x3  (aj,ajinv,det) 
c 

c  ***  compute  global  derivatives  and  strain-nodal  displacement  matrix 
c 

do  160  i=l,n 

derivg (1, i) =ajinv (1, 1) *derivl (1, i) +ajinv(l, 2) *derivl (2, i) 
derivg { 2 , i ) =a j inv (2,1) *derivl ( 1 , i ) +a j inv (2,2) *derivl ( 2 , i ) 
der i vg ( 3 , i ) =a j inv (3,1)*  derivl ( 1 , i ) +a j inv (3,2) *deri vl ( 2 , i ) 
derilg(l,i)=ajinv(l,3) *hzdt 
derilg (2, i) =ajinv (2, 3) *hzdt 
derilg(3,i)=ajinv(3,3) *hzdt 
160  continue 

c 

do  170  i=l, 6 
do  170  j-1,24 

bmtx (i, j )=0 . 0 

170  continue 

c 

do  180  i=l,n 

il= (i-1) *6+1 

i2=il+l 

i3=i2+l 

i4=i3+l 

i5=i4+l 

i6=i5+l 

gkl=derivg (1, i) *hz+shapef (i) * derilg (l,i) 
gk2=derivg (2, i) *hz+shapef (i) * derilg (2, i) 
gk3=derivg (3, i) *hz+shapef (i) * derilg (3>i) 
c 

bmtx ( 1 , il ) =derivg ( 1 ,  i ) 
bmtx(l,i4)=gkl* (-v2 (1) ) 
bmtx ( 1 , i5 ) =gkl* vl ( 1 ) 
bmtx(l,i6)=gkl*v3 (1) 
bmtx ( 2 , i2 ) =derivg ( 2 ,  i ) 
bmtx(2,i4)=gk2*(-v2{2) ) 
bmtx (2, i5) =gk2*vl (2) 
bmtx (2, i6) =gk2*v3 (2) 
bmtx ( 3 , i3 ) =derivg ( 3 , i ) 
bmtx(3,i4)=gk3* (-v2 (3) ) 
bmtx (3, i5) =gk3*vl (3) 
bmtx (3, i6) =gk3*v3 (3) 
bmtx ( 4 , il ) =derivg ( 2 ,  i ) 
bmtx ( 4 , i2 ) =derivg ( 1 ,  i ) 
bmtx (4, i4) =gk2* (-v2 (1) ) +gkl* (-v2(2) ) 
bmtx (4, i5) =gk2*vl (1) +gkl*vl (2) 
bmtx(4,i6)=gk2*v3 (1) +gkl*v3 (2) 
bmtx ( 5 , i2 ) =derivg (3,  i) 
bmtx ( 5 , i3 ) =derivg ( 2 , i ) 
bmtx (5, i4) =gk3* (-v2 (2) ) +gk2* (-v2(3) ) 
bmtx (5, i5) =gk3*vl (2) +gk2*vl (3) 
bmtx ( 5 , i 6 ) =gk3* v3 ( 2 ) +gk2  * v3 ( 3 ) 
bmtx (6, il) =derivg (3, i) 
bmtx ( 6, i3) =derivg (1, i) 
bmtx (6, 14) =gk3* (-v2 (1) ) +gkl* (-v2(3) ) 
bmtx (6, i5) =gk3*vl (1) +gkl*vl (3) 
bmtx (6, 16 ) =gk3*v3 (1) +gkl*v3 (3) 
continue 


180 
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c 

c  ***  calculate  transpose  of  bmtx 
c 

do  190  i=l,6 
do  190  j=l,neldof 

190  bmtxt (j,i)=bmtx(i, j) 

c 

weight  =  wx  *  wy  *  wz 

detwt  =  det  *  weight 
icount^icount+l 
tmp (icount ) =detwt 
c 

c  ***  calculate  strain 
c 

call  mmult (6,24, l,bmtx, edisp, estrainp, 0,1.0) 
c 

c  ***  Transform  strain  to  local  coordinates  system 
c 

call  mmult (6, 6, l,rot6, estrainp, estrain, 0,1.0) 
c 

c  ***  Subtract  any  prior  plastic  strain 
c 

do  200  i=l,6 

200  estrain(i)  =  estrain(i)  -  stssh4 (12+i, issts) 

c 

stressO  =  stssh4 (19, issts) 

Phi  =  stssh4 (20, issts) 
plastrn  =  stssh4 (21, issts) 
c 

c  ***  Calculate  stress  using  plane-strain  formulas 
c 

estress (1)  =  elast/ (1 . 0-poiss**2) * (estrain (1)  +  poiss*estrain (2) ) 
estress (2 )  =  elast/ (1 . 0-poiss**2) * (poiss*estrain (1)  +  estrain(2)) 
estress(3)  =  elast  *  estrain (3) 
estress (4)  =  shear  *  estrain (4) 
estress  (5)  =  5.0  /6.0  *  shear  *  estrain(5) 
estress (6)  =  5.0  /6.0  *  shear  *  estrain (6) 
c 

c  ***  Elastic-Plastic  Calculations 
c 

p  =  - (estress (1)  +  estress (2)  +  estress (3))  /  3.0 
c 

do  210  i=l, 6 

210  s(i)  =  estress(i)  +  p  *  delt(i) 

q  =  sqrt (1.5  *(s(l)**2  +  s(2)**2  +  s(3)**2  +  2.0  *  s(4)**2  + 
1  2.0  *  s  (5) **2  +  2.0  *  s (6) **2) ) 

F  =  (q  /  stressO) **2+2. 0*ql*Phi*cosh(-1.5*q2*p/stress0)- 
1  (1.0+q3*  Phi**2) 

c 

if  (F.gt.0)  then 
iter  =  0 
dep  =  0.0 
deq  =  0.0 
c 

c  **  Begin  Iteration  Loop  ** 
c 

cl  =  3.0*q2/(2.0*stress0**2) 

Fql  =  2.0  *  q  /  (stress0**2) 

Fq2  =  2.0  /  (stress0**2) 
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Fpl  =  2.0  *  ql  *  Phi  *  -cl  *  sinh{-cl  *  p) 

Fp2  =  2.0  *  ql  *  Phi  *  cl**2  *  cosh (-cl  *  p) 

Fsl  =  -q  ** 2  /  (2.0  *  stress0**3)  +  2 . 0*ql*Phi*cl* 
p*sinh(-cl*p) 

Fsp  -  2.0*ql*Phi* (~cl**2  *  p*cosh (-cl*p)  + 

1  cl/stressO  *  sinh(-cl*p)) 

Fsq  =  -q  /  stress0**3 

stressl  =  Syp 

esl  =  stressO/elast 

el  =  -p*dep/ ( (1-Phi) *stressO) 

psl  =  plastrn  +  el  +  esl 

iflag  =  0 

do  230  i=l,iseg 

if  (iflag  .eq.  0)  then 
if (psl  .It.  strain (i+1 ) )  then 
stressl  =  stressl  +  slope (i) * (psl-strain (i) ) 
isegl  =  i 
iflag  =1 
else 

stressl  =  stressl  +  slope (i) * (strain (i+1) -strain (i) ) 
endif 
endif 

230  continue 
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c 

c  ***  Numerical  Error  Trap 
c 

if  ((stressl  .eq.  stressO)  .or.  (abs(el)  .eq.  0.0)) 

then 

dsdep  =  slope (isegl) 

else 

dsdep  =  (stressl-stressO) /el 

endif 

if  (abs (dsdep)  .gt.  elast)  then 
dsdep  =  slope (isegl) 

endif 

c 

stressl  -  Syp 

el  =  q  *  deq  /  ((1-Phi)  *  stressO) 
psl  -  plastrn  +  el 
iflag  =  0 
do  240  i=l,iseg 
if  (iflag  .eq.  0)  then 
if (psl  .It.  strain (i+1))  then 
stressl  =  stressl  +  slope (i) * (psl-strain (i) ) 
isegl  =  i 
iflag  «  1 
else 

stressl  =  stressl  +  slope (i) * (strain (i+1) - 

strain(i)) 

endif 

endif 

240  continue 

c 

c  ***  Numerical  Error  Trap 
c 

if  ((stressl  .eq.  stressO) .or . (abs (el)  .eq.  0.0))  then 
dsdeq  =  slope (isegl) 

else 

dsdeq  -  (stressl-stressO) /el 

endif 

if  (abs (dsdeq)  .gt.  elast)  then 
dsdeq  =  slope (isegl) 

endif 

c 

c  ***  Using  only  stressO  derivatives 
c 

A  (1,1)  -  bulk  *  Fpl+Fsl*dsdep 
A(l,  2)  =  -3. 0*shear*Fql+Fsl*dsdeq 

A(2, 1)  =  Fql+deq* (bulk*Fp2+Fsp*dsdep) +dep*Fsq*dsdep 
A  (2, 2)  =  Fpl+dep* (-3 . 0*shear*Fq2+Fsq*dsdeq) +deq*Fsp*dsdeq 
c 

B ( 1 )  =  -1.0  *  F 

B (2)  =  -1.0  *  dep  *  Fql  -  deq  *  Fpl 
c 

call  inv2x2 (A,Ainv,d) 
c 

call  mmult (2,2, l,Ainv,B, corr, 0,1.0) 
c 

Phi  =  stssh4 (20, issts) 
plastrn  =  stssh4 (21, issts) 
stressO  =  stssh4 (19, issts) 
c 

dep  =  dep  +  corr ( 1 ) 
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deq  =  deq  +  corr(2) 
c 

degrowth  =0.0 
do  250  i=l,5 

ind  =  iplane(i) 

dstrnplst (ind)  =  dep  *  delt(ind)  /  3.0  + 

1  deq  *-3.0  *  s (ind)  /  (2.0  *  q) 

250  continue 

c 

dstrnplst (3)  =  dep  /  3.0  +  deq  *  estress(3)  /  q 
c 

do  260  i=l, 3 

260  '  degrowth  =  degrowth  +  dstrnplst (i) 

c 

c  **  Calculate  New  stress,  p,  q,  and  F 
c 

estress(l)  =  elast/ (1 . 0-poiss**2) * (estrain (1)  -  dstrnplst(l)  + 

1  poiss* (estrain (2) -dstrnplst (2) ) ) 

estress(2)  =  elast/ (1 . 0-poiss**2) * 

1  (poiss* (estrain (1) -dstrnplst (1) ) testrain (2) -dstrnplst  (2) ) 

estress(3)  =  elast  *  (estrain (3) -dstrnplst (3) ) 
estress(4)  =  shear  *  (estrain (4) -dstrnplst (4) ) 
estress{5)  =  5.0  /6.0  *  shear  *  (estrain (5) -dstrnplst (5) ) 
estress(6)  =  5.0  /6.0  *  shear  *  (estrain ( 6) -dstrnplst (6) ) 
p  =  (estress(l)  +  estress(2)  +  estress(3))  /  -3.0 
c 

do  270  i=l, 6 

270  s(i)  =  estress(i)  +  p  *  delt(i) 

c 

deltep  =  (q*deq  -  p  *  dep) / ((1-Phi) *stress0) 
c 

c  Prevent  reduction  in  plastic  strain 

if  (deltep  .It.  0.0)  then 
deltep  =0.0 

endif 

c 

plastrn  =  plastrn  +  deltep 
if  (fn  .ne.  0.0)  then 

anucl  =  fn/ (sn*2. 50662827463) *exp (-.5* ( (plastrn-en) /sn) **2) 
else 

anucl  =  0.0 
endif 
c 

Phil  =  Phi  +  (1  -  Phi)  *  de growth  +  anucl  *  deltep 
c 

c  Phi  cannot  decrease 

if  (Phil  .gt.  Phi)  Phi  =  Phil 
c 

c  Turn  off  Void  Effects  if  ql  =  0.0 

if  (ql  .eq.  0.0)  Phi  =  0.0 
c 

c  Calculate  Strain  Hardening 

elastrn  =  stressO/elast 
stressO  =  Syp 
iflag  =  0 
do  280  i=l,iseg 

if  (iflag  .eq.  0)  then 
psl  =  plastrn  +  elastrn 
if(psl  .It.  strain (i+1) )  then 
stressO  =  stressO  +  slope (i) * (psl-strain (i) ) 
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iflag  =  1 
else 

stressO  =  stressO  +  slope (i) * (strain (i+1) -strain (i) ) 
endif 
endif 

280  continue 

c 

q  =  sqrt (1.5  *(s(l)**2  +  s(2)**2  +  s(3)**2  +  2.0  *  s(4)**2  + 
1  2.0  *  s(5)**2  +  2.0  *  s  (6) **2) ) 

c 

F  =  (q  /  stressO) **2+2. 0*ql*Phi*cosh (-1 . 5*q2*p/stress0) - 
1  ( 1 . 0+q3*  Phi* *2 ) 

c 

iter  =  iter  +  1 
c 

if  (iter. gt. 500)  then 

write (*,*)  'Failed  to  Converge1 
stop 

endif 

if (abs (F) .gt.eps)  goto  220 
c 

endif 

c  ***  End  of  Plastic  Deformation  Section 

c 

c 

c  ***  Transform  stress  back  to  global  coordinates 
c 

call  mmult (6,6,1, rot6t, estress, estressg, 0,1.0) 
c 

c  ***  find  eforce  for  this  gauss  point  and  sum  to  element  force 
c 

call  mmult (24, 6, l,bmtxt, estressg, eforce,  1, detwt ) 
c 

c  ***  save  stress  and  strain 
c 

do  290  i  =  1,6 

stssh4 (i, issts) =estress ( i ) 

stssh4 (i+6,issts)=estrain(i)  -  dstrnplst(i) 
stssh4 (i+12, issts) =stssh4 (i+12, issts)  +  dstrnplst (i) 
290  continue 

c 

stssh4 (19, issts)  =  stressO 
stssh4  (20,  issts)  ==  Phi 
stssh4 (21, issts)  -  plastrn 
c 

c  ***  end  of  integration  loop 
c 

300  continue 

c 

c  ***  Hourglass  Control  Correction 
c 

do  310  i  =  1,24 

310  eforce (i)  -  eforce (i)  +  hour  *  (eforceh(i)  -  eforce (i) ) 

c 

c  ***  Transform  Moments  back  to  Global  Coordinate  system 
c 

do  320  i=l, 4 
do  320  j=l,3 

eforceg ( (i-1) *6+j ) =eforce ( (i-1) *6+j ) 
eforceg ( (i-1) *6+3+j ) =eforce ( (i-1) *6+4 ) *tl ( j , 1) 
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c 

c 

c 


1 
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continue 


+eforce ( (i-1) *6+5) *tl  ( j , 2) 
+eforce ( (i-1) *6+6) *tl ( j , 3) 


320 

*** 


340 

c 

c 

350 

c 

q  ★  -k 

c 


360 

c 

c 

370 

c 


compute  lumped  mass  matrix 


sum=0 . 0 

do  340  i=l,  ngauss 
sum=sum+tmp  ( i ) 

ems=sum/n 

do  350  i=l,  neldof 
emass (i) =ems*densty 

assemble  mass  and  force  matrices  into  system  matrix 
do  360  i-1,  neldof 

sforce (lindex (i) )  =  sforce (lindex (i) )  +  eforceg(i) 
smass (lindex (i) )  ~  smass (lindex (i) )  +  emass (i) 
continue 

end  of  element  loop 


continue 


return 

end 
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